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Abstract 

The  age  of  permafrost  is  closely  linked  to  the  time  required  for  soil  systems  to 
freeze,  since  the  permafrost  must  be  at  least  as  old  as  the  formation  time.  Cycles 
of  freeze-thaw  will  complicate  the  relation  between  the  freeze  rate  and  the  age. 
A  model  based  on  pure  conduction  heat  transfer  with  freeze-thaw  is  used  to 
predict  the  time  required  for  a  given  thickness  of  permafrost  to  develop,  either 
heterogenetically  or  syngenetically.  The  formation  time  is  a  function  of  the  long¬ 
term  geothermal  gradient  (initial  temperature  of  the  thawed  soil),  the  ratios  of 
the  frozen  to  thawed  thermal  properties,  and  the  temperature  history  of  the  upper 
surface  of  the  permafrost  (higher  than  the  air  temperature).  The  simple  theory 

allows  universal  graphs  to  be  produced  that  predictthe  formation  time  for  a  given 

thickness  of  permafrost.  Realistic  soil  property  ratios  and  paleotemperature 
scenarios  will  then  lead  to  estimates  of  the  formation  time  of  permafrost  for  a 
specific  site.  The  model  indicates  that  deep  permafrost  (more  than  1500  m) 
requires  formation  times  on  the  order  of  the  complete  Quaternary  Period. 
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PREFACE 
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Development  Program,  Deep  Permafrost  Borehole  Sites  in  Alaska,  2  6.00  78  203,  and  DA  Project 
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NOMENCLATURE 


p  8/X 


c  specific  heat 

C12  Cl  /C2 

C  pc,  volumetric  specific  heat 

/  fraction  of  geothermal  energy  used  for 
melt 

a2i(l  +  cj) 

P[(P  +  2)(J  +  2] 

g'  dg/da 

G  geothermal  gradient  (slope  of  initial  tem¬ 
perature  distribution) 

k  thermal  conductivity 
kg  thermal  conductivity  of  soil  solids 
k[2  k[lk2 

I  latent  heat  of  solidification 
L  p£,  volumetric  latent  heat 
M  nfiuhMtc 
qg  geothermal  energy  flow 
2s»  Ql  sensible  and  latent  heats 
S  soil  saturation  parameter 

Stefan  Number 

t  time 

k  characteristic  time  for  perma¬ 

frost  temperature  changes 

time  to  complete  melt  of  permafrost 
T  temperature 

U  deposition  rate  at  surface  of  permafrost 
X  depth  coordinate 
Xvv  volumetric  water  content 
X  permafrost  thickness 
a  thermal  diffusivity 


Subscripts 

f  solidification  or  freeze  value 
e  equilibrium  value 
i  initial  value 
o  initial  surface  value 
s  surface 
u  thawed 

1 ,2  frozen  and  thawed,  respectively 
oo  steady-state  value 


Permafrost  Formation  Time 


VIRGIL  J.  LUNARDINI 


INTRODUCTION 

The  age  of  permafrost  is  of  interest  to  biologists,  geophysicists  and  engineers.  Clearly,  permafrost  must 
be  at  least  as  old  as  the  time  it  took  for  it  to  form;  thus,  the  formation  time  of  permafrost  can  be  considered 
its  minimum  age.  A  volume  of  permafrost  can  be  much  older  than  this  since  it  may  exist  for  many  years  af¬ 
ter  formation.  This  report  will  examine  the  formation  time  of  permafrost  using  a  pure  heat  conduction  mod¬ 
el.  As  we  shall  see,  the  surface  temperature  history  of  the  soil  mass  is  critical  for  any  prediction  of  the  per¬ 
mafrost  formation  time.  Since  the  formation  time  of  permafrost  is  expected  to  be  on  the  order  of  millennia, 
it  is  necessary  to  examine  the  geophysical  record  to  obtain  some  bounds  on  realistic  surface  temperatures 
that  the  Earth  has  experienced  during  the  time  when  permafrost  was  growing.  First,  we  will  discuss  perma¬ 
frost  and  paleotemperature  scenarios,  then  we  will  formulate  a  mathematical  model  of  permafrost  growth, 
and,  finally,  we  will  examine  some  predictions  by  the  model  of  permafrost  formation  times. 

Permafrost  is  a  widespread  phenomenon  that  has  been  and  still  is  greatly  misunderstood.  The  term 
“permafrost”  is  generally  attributed  to  S.W.  Muller  (1945),  who  apparently  coined  the  name  in  place  of  the 
more  awkward  terms:  permanently  frozen  ground  or  permanent  frost.  Bryan  (1946)  suggested  the  term 
“pergelisol,”  but  this  has  not  been  adopted  except  in  the  French  literature.  In  order  to  understand  the  con¬ 
cept,  let  us  look  at  a  general  definition  given  in  Lunardini  (1981a): 

Permafrost  describes  the  thermal  condition  of  earth  materials  (sand,  glacial  till,  organic  matter,  etc.) 
when  their  temperature  remains  at  or  below  32°F  (0°C)  continuously  for  a  significantly  long  time, 
but  not  necessarily  for  an  entire  geological  period.  It  does  not  include  earth  materials  that  drop  be¬ 
low  32°F  during  one  winter  and  remain  below  32°F  through  the  following  summer  and  into  the  next 
winter,  although  for  practical  engineering  purposes  such  materials  may  be  included. 

Clearly,  permafrost  is  not  so  much  a  material  as  it  is  the  thermal  state  of  ordinary  soil  systems.  It  does 
not  include  systems  that  are  at  or  below  0°C,  but  contain  no  earth  materials,  e.g.,  ice  caps,  glaciers  and  ice¬ 
bergs.  There  is  no  agreement  on  the  minimum  time  during  which  the  material  must  remain  below  0°C  to 
qualify  as  permafrost.  Soils  that  freeze  during  an  exceptionally  severe  winter  and  survive  for  1  or  2  more 
years  are  called  “pereletoks,”  and  often  are  not  classified  as  permafrost  (Swinzow  1969). 

The  existence  of  permafrost  is  a  result  of  the  history  and  the  present  state  of  the  energy  balance  at  the 
Earth’s  surface — measured  by  the  surface  temperature — and  the  deep  Earth  heat  flow.  If  permafrost  exists 
and  the  net  yearly  gain  of  energy  by  the  entire  permafrost  volume  is  equal  to  the  net  loss  of  energy,  then  the 
permafrost  will  remain  stationary,  while  an  excess  heat  gain  over  heat  loss  will  result  in  a  net  loss  of  frozen 
material.  Given  the  same  energy  balances,  however,  i.e.,  net  gain  of  energy  over  the  year,  one  region  may 
have  permafrost  (albeit  degrading)  while  another  will  not.  This  is  ascribable  to  the  thermal  history  of  the 
frozen  ground  in  the  two  areas.  Though  both  are  losing  or  have  lost  permafrost,  one  region  may  have  start¬ 
ed  with  a  larger  volume  of  permafrost  than  the  other.  Thus,  the  present  energy  flow  conditions  may  be  such 
that  permafrost  cannot  exist  in  one  region,  whereas  it  will  subsist  in  another  area,  although  in  a  receding 
form  often  referred  to  as  “relic  permafrost.” 


In  this  sense  previous  glaciation  has  very  likely  played  an  important  role  in  the  present  existence  of  per¬ 
mafrost  in  marginal  areas.  It  is  safe  to  conclude  that  little,  if  any,  permafrost  exists  beneath  nonpolar  gla¬ 
ciers,  but  once  they  withdraw,  permafrost  may  rapidly  form  and  grow.  So,  previously  glaciated  regions  will 
show  a  lesser  volume  of  frozen  ground  than  unglaciated  regions  with  similar  climatic  histories.  In  this  regard 
it  is  significant  that  Canada  was  heavily  glaciated  while  Russia  (Siberia)  had  little  permanent  glaciation. 
Thus,  the  permafrost  thickness  in  Siberia  is  much  greater  than  in  Canada,  although  the  climates  are  similar. 

Since  the  energy  balance  involves  meteorological  conditions,  surface  vegetation,  topography  and  soil 
conditions,  we  can  anticipate  that  we  will  find  no  simple  correlation  for  the  existence  of  permafrost  in  the 
marginal  or  discontinuous  regions.  Permafrost  does  exist  far  south  of  the  usually  accepted  limits  in  scattered 
patches  almost  always  associated  with  high  altitudes  and,  thus,  microclimates  similar  to  the  usual  permafrost 
regions. 

Origin  and  existence  of  permafrost 

In  discussing  the  present  distribution  of  permafrost,  questions  often  arise  concerning  its  origin  and  age. 
These  two  concepts  should  be  clearly  differentiated,  since  they  deal  with  two  separate  phenomena.  The  age 
of  a  particular  deposit  of  permafrost  is  the  time  that  has  elapsed  since  the  freezing  of  the  soil  system.  Actual¬ 
ly,  it  may  be  very  difficult  or  impossible  to  determine  this  age  because  thawing  and  freezing  may  cycle  at 
long  intervals  and  different  frequencies  in  different  regions  of  the  Earth.  Thus,  the  ages  of  two  “similar”  de¬ 
posits  of  permafrost  may  be  quite  different.  In  this  regard,  the  presence  of  preserved  animal  remains  may  be 
a  reliable  clue  to  the  age  of  a  deposit  of  frozen  ground.  The  age  of  permafrost  is  a  question  of  significance 
and  may  be  useful  to  paleontologists,  paleobotanists,  etc.  The  present  thermal  state  of  the  permafrost— tem¬ 
perature,  degradation,  aggradation,  etc. — is  of  interest  to  engineers. 

The  origin  of  permafrost  involves  the  question  of  the  conditions  under  which  it  can  form  and  grow.  These 
same  conditions  will  explain  the  present  existence  of  permafrost  at  a  given  location.  As  the  conditions  for 
the  origin  of  permafrost  are  dynamic,  it  is  certainly  possible  that  areas  now  lacking  permafrost  once  had 
these  underlying  frozen  strata,  and  that  the  present  regions  of  permafrost  could  once  have  been  thawed.  In 
other  words,  the  present  existence  of  permafrost  depends  upon  two  things:  the  proper  energy  exchange  con¬ 
ditions  and  the  thermodynamic  state  of  the  permafrost  mass  itself.  The  first  of  these  conditions  has  little  to 
do  with  past  climatic  conditions,  but  the  second  is  a  function  of  the  complete  thermal  history  of  the  perma¬ 
frost  and  is  thus  related  to  past  climate.  In  this  sense,  it  is  incorrect  to  describe  permafrost  simply  as  a  legacy 
of  the  last  great  ice  age.  It  is  possible  that  some,  perhaps  most,  permafrost  had  its  origin  at  the  beginning  of 
the  Pleistocene  era  (Brown  1964),  but  this  should  not  imply  that  the  intervening  thermal  conditions  were 
without  significance. 

From  the  above  discussion,  it  should  be  clear  that  the  formation  and  existence  of  permafrost  are  related  to 
the  present  and  historical  conditions  of  energy  exchange  between  the  soil  and  the  atmosphere.  Nevertheless, 
it  is  not  possible  to  state  these  conditions  in  a  simple,  precise  manner  that  will  allow  us  to  define  unique  per¬ 
mafrost  indicators. 

Over  a  sufficiently  long  time  span,  the  energy  exchange  will  be  periodic,  and,  averaged  over  a  number  of 
periods,  the  net  energy  flow  for  a  given  soil  volume  will  determine  its  thermodynamic  state.  Dynamic  equi¬ 
librium  of  the  energy  flows  may  exist  such  that  the  soil  is  perennially  frozen  or  the  thermodynamic  state  of 
the  soil  may  be  varying  because  of  an  imbalance  of  the  cyclic  energy  flows.  The  original  formation  of  per¬ 
mafrost  depends  upon  a  net  periodic  (yearly)  loss  of  energy  from  the  soil  volume  that  must  persist  for  many, 
many  years  for  the  permafrost  to  attain  great  thicknesses.  The  maintenance  of  the  present  thermodynamic 
state  of  permafrost  requires  only  that  the  net  energy  flow  averaged  over  a  number  of  years  be  zero. 

Paleotemperatures 

The  thermal  history  of  permafrost  is  greatly  influenced  by  the  long-term  temperature  variations  experi¬ 
enced  at  its  upper  surface.  The  relative  mean  global  temperature  deduced  from  oxygen  isotope  data  is  shown 
for  the  past  180  million  years  in  Figure  1  (Eddy  and  Bradley  1991).  There  was  probably  little  or  no  perma¬ 
frost  prior  to  the  late  Tertiary  Period  and  certainly  none  for  100  million  years  prior  to  the  long-term  cooling 
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Figure  L  Paleotemperature  varia¬ 
tions  (after  Eddy  and  Bradley  1991). 
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Figure  2.  Pleistocene  temperature  fluc¬ 
tuations  (after  Folland  et  aL  1990). 


that  began  35  million  years  before  present.  During  the  Pliocene,  some  2-5  million  years  before  present,  the 
temperature  oscillated  0-4°C  above  the  present  values.  Permafrost  was  probably  present  at  locations  that 
now  have  mean  temperatures  less  than  -14°C  but  with  greatly  reduced  thickness  and  quite  variable  temporal 
existence. 

The  Quaternary  Period,  which  includes  the  Pleistocene  (about  1  million  years)  and  the  Holocene  (present 
to  about  10,000  years),  remains  a  time  of  greatly  reduced  temperatures  and  massive  glaciations.  The  period 
is  marked  by  ice  ages  of  100,000  to  120,000  years  duration,  interrupted  by  interglacials,  and  with  increas¬ 
ingly  severe  minimum  temperatures  and  temperature  drops.  We  are  now  in  the  Holocene  interglacial  (which 
is  actually  cooler  than  the  previous  interglacial,  the  Sangamon  or  Eem)  and  the  record  would  seem  to  indi¬ 
cate  the  next  significant  temperature  move  should  be  downward  with  a  new  ice  age  evolving.  Figure  2 
shows  the  temperature  record  (again  from  isotope  analysis)  for  the  past  million  years  or  so  (Folland  et  al. 
1990).  We  note  graphically  the  record  of  at  least  four  or  five  glacial  periods,  with  the  last  one  ending  some 
12,000  years  ago  in  North  America — the  Wisconsin.  The  temperatures  variations  from  present  values  swung 
from  highs  of  +3°C  during  interglacials  to  lows  of  about  -10°C  during  the  glacial  maximums  of  the  last  mil¬ 
lion  years.  The  overall  trends  reveal  very  rapid  temperature  rises  over  time  spans  of  approximately  12,000- 
13,000  years,  followed  by  less  rapid  temperature  drops  over  20,000  years,  followed  by  a  period  of  about 
100,000  years  of  gradually  decreasing  temperatures  with  interglacial  temperature  rises  of  4®C. 

The  temperature  history  of  the  past  160,000  years  has  been  quantified  using  the  deep  glacial  ice  cores 
taken  in  Greenland  and  Antarctica.  Figure  3  shows  the  temperature  variations  from  present  values,  for  the 
time  of  the  last  great  ice  age,  taken  from  the  isotope  analysis  of  the  Vostok  (Antarctic)  core  (Jouzel  et  al. 
1987).  The  maximum  temperature  excess  of  the  past  160,000  years  was  about  3°C  and  thus  only  the  discon¬ 
tinuous  permafrost  zones  would  have  been  in  danger  of  disappearing,  although  there  was  continual  change 
in  the  permafrost  thickness.  If  we  assume  that  the  temperature  variations  obtained  from  the  ocean  isotope 
and  ice  core  analyses  are  global,  then  the  paleotemperature  history  at  a  specific  site  can  be  reconstructed  us¬ 
ing  the  present  temperature  data.  This  has  been  done  for  Prudhoe  Bay  by  Osterkamp  and  Gosink  (1991)  as 
shown  in  Figure  4  and  for  Fast  Siberia  by  Maximova  and  Romanovsky  (1988)  (Fig.  5).  Sun  and  Li  (1988) 
also  presented  a  quantitative  model  of  temperature  fluctuations  during  the  last  ice  age  in  northern  China 
(Fig.  6). 

The  isotopic  data  all  indicate  that  the  temperature  at  a  given  site  could  have  dropped  as  much  as  10  to 
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Temperature  Difference  From  Present  (°C) 


Figure  3.  Vostok  (Antarctica)  ice  core  tem¬ 
perature  inferences  (after  Jouzel  et  aL  1987) 


c.  Based  on  Brigham  and  Miller  ( 1983)  data 
for  Barrow,  Alaska. 


Figure  4.  Permafrost  surface  paleotemperature  model  for  Prudhoe  Bay,  Alaska  (after  Osterkamp 
and  Go  sink  1991). 
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Figure  5.  Permafrost  surface  paleotem- 
perature  model  for  East  Siberian  location 
(after  Maximova  and  Romanov  sky  1988). 


Figure  6.  Reconstructed  temperature 
for  northern  China  during  last  ice  age 
(after  Sun  and  Li  1988). 
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a.  Based  on  Vostok  (Antarctica)  ice  core  data,  b.  Based  on  Brigham  and  Miller  (1983)  data, 
AT  =  -5. 29  AT  =  -2. 54  °C. 


Figure  7,  Mean  paleotemperature  departure,  Prudhoe  Bay,  Alaska,  during  one  glacial  cycle. 


12°C  below  the  present  values  for  varying  periods  of  time.  While  extreme  temperature  drops  at  some  sites 
may  have  been  significantly  greater  than  these  values,  there  is  no  convincing  evidence  of  this.  Temperature 
variations  of  8  to  12°C  during  the  last  glaciation  have  also  been  reported  from  the  Greenland  ice  cores 
(Dansgaard  and  Oeschger  1989).  Folland  et  al.  (1990)  note  that  global  temperatures  underwent  5-7°C  vari¬ 
ations,  with  changes  as  great  as  10-15°C  at  middle  and  high  latitude  regions  of  the  Northern  Hemisphere. 
During  the  Eemian  interglacial,  temperatures  in  Siberia,  Canada  and  Greenland  may  have  increased  by  4- 
8°C.  Thus,  it  would  seem  prudent  to  use  extremal  paleotemperature  excursions  on  this  order  of  magnitude 
(10-12°C)  for  the  development  of  permafrost  formation  models.  We  cannot  ascribe  the  rapid  growth  of 
deep  permafrost  to  extraordinarily  low  temperatures  that  are  beyond  the  ranges  we  have  mentioned  above. 
Figure  7  shows  two  examples  of  average  paleo temperatures  for  a  glacial  cycle  made  up  of  glacial  and  inter- 
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Figure  8.  Paleotemperature  history  at  Macken¬ 
zie  Delta,  Canada  (after  Allen  et  al  1988). 
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glacial  intervals.  Figure  7a  is  based  on  the  Vostok  ice  core  while  Figure  7b  is  calculated  from  information 
given  by  Brigham  and  Miller  (1983)  for  Prudhoe  Bay,  Alaska. 

There  have  been  surprisingly  few  systematic  thermal  studies  of  the  origin  of  permafrost  and  the  total 
time  required  for  its  formation.  Osterkamp  and  Gosink  (1991)  studied  the  response  of  permafrost  thickness 
to  surface  temperature  variations.  They  used  quasi-steady  and  numerical  models  to  predict  the  position  of 
the  permafrost  bottom,  with  arbitrary  initial  permafrost  thicknesses,  but  were  interested  in  the  inverse  prob¬ 
lem  of  deducing  paleotemperatures  from  the  present  permafrost  data  at  Prudhoe  Bay,  Alaska.  Allen  et  al. 
(1988)  used  a  quasi-steady  model  (which  will  inherently  underpredict  the  growth  time)  for  the  same  pur¬ 
pose  in  the  Mackenzie  Delta  region  of  Canada,  using  still  another  paleotemperature  history  (Fig.  8).  They 
assumed  that  Illinois  permafrost  was  completely  melted  during  the  Sangamon  interglacial,  although  this  is 
highly  unlikely  given  the  time  available,  unless  Illinoisian  glaciation  limited  the  permafrost  thickness  to 
modest  values.  Nevertheless,  they  used  an  initial  permafrost  thickness  and  predicted  about  800  m  of  growth 
in  75,000  years.  Their  thermal  model  is  based  on  the  work  of  Lachenbruch  et  al.  (1982)  dealing  with  the 
rate  of  thaw  of  thick  permafrost  zones  during  several  millennia.  Romanovsky  et  al.  (1988)  noted  qualitative 
aspects  of  the  origin  and  disappearance  of  permafrost  in  the  Transbaikal  area  of  Russia.  Katasonov  (1988) 
used  “cryogenic  structures”  to  argue  for  the  origin  of  permafrost  early  in  the  Quaternary  and  its  persistence 
up  to  the  present. 

In  this  report  an  attempt  is  made  to  predict  the  rate  of  permafrost  formation,  starting  with  no  permafrost, 
i.e.,  its  origin,  using  a  simple  conduction  model.  If  the  soil  forming  the  permafrost  exists  before  freezing 
starts,  the  growth  is  heterogenetic.  When  the  permafrost  forms  as  the  soil  material  is  gradually  deposited  at 
the  surface,  the  permafrost  is  said  to  have  a  syngenetic  origin.  The  thermal  conditions  for  each  type  of 
growth  will  be  examined. 

THEORY 

The  solution  to  conductive  heat  transfer  problems,  with  solidification  phase  change,  has  interested  engi¬ 
neers  and  mathematicians  for  over  a  century.  These  problems  (often  referred  to  as  Stefan  problems)  are  in¬ 
herently  nonlinear  and  solution  methods  are  very  restricted.  A  classical  solution  for  the  case  of  a  constant 
temperature,  semi-infinite  medium  that  undergoes  a  step  change  of  surface  temperature  was  given  by  Neu¬ 
mann  (1860)  and  expanded  upon  by  Carslaw  and  Jaeger  (1959);  it  is  called  the  Neumann  solution.  Tao 
(1978)  extended  the  similarity  technique  of  Neumann  to  the  semi-infinite  slab  with  arbitrary  initial  temper¬ 
ature.  This  is  precisely  the  solution  we  need,  but  unfortunately  this  exact  solution  is  such  that  numerical 
computations  are  extremely  difficult  because  of  transient  functions  that  require  an  increasing  number  of  se¬ 
ries  terms  as  time  increases.  For  permafrost  formation  the  time  scales  are  so  huge  that  Tao’s  solution  is  im¬ 
practical  and  cannot  be  used.  Like  the  exact  solution  of  Lozano  and  Reemsten  (1981),  for  flux  boundary 
conditions,  Tao’s  solution  is  perhaps  best  used  numerically  to  verify  the  accuracy  of  approximate  and  nu¬ 
merical  solutions  or  for  short-time  solutions. 

The  search  for  practical  solutions  for  engineering  design  has  led  to  some  convenient  approximate  solu¬ 
tion  methods  for  Stefan  problems.  The  heat  balance  integral  technique  solves  the  energy  equation  on  aver- 
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age  over  a  space  volume,  instead  of  at  each  point  of  space  (Goodman  1958).  The  concept  is  identical  to  the 
well-known  momentum  integral  technique  of  fluid  mechanics,  sometimes  referred  as  the  Karman-Pohlhaus- 
en  method  (Schlichting  1968).  The  method  has  been  applied  successfully  to  constant  initial  temperature 
problems  of  the  semi-infinite  slab  (Lunardini  and  Varotta  1981)  as  well  as  the  cylindrical  geometry  (Lunar- 
dini  1980).  A  modification  of  the  integral  method  utilizing  a  single  integration  over  an  entire  nonconstant 
property  volume  has  yielded  accurate  solutions  (Yuen  1980,  Lunardini  1981b, 1982). 

The  integral  solution  has  been  used  for  a  problem  with  variable  initial  temperature  distributions,  but  the 
results  were  limited  to  shallow  freeze  depths  (Lunardini  1984).  This  report  presents  an  approximate  solution 
to  the  modified  Neumann  problem  for  which  a  linear  initial  temperature  distribution  exists.  Such  an  initial 
temperature  is  common  for  soil  systems  with  a  geothermal  temperature  gradient  and  is  directly  applicable  to 
the  question  of  permafrost  formation  rates. 

Heterogenetic  freeze  relations 

Figure  9  shows  the  case  of  an  infinite  layer  of  soil  with  a  linear  initial  temperature  distribution  (G  repre¬ 
sents  a  geothermal  gradient).  The  soil  is  assumed  to  be  homogeneous  and  conduction  is  the  only  mode  of 
energy  transfer.  At  zero  time  the  surface  temperature  drops  to  and  is  held  constant  while  freezing  com¬ 
mences.  At  any  time  t  >  0,  there  is  a  frozen  layer,  called  layer  1  (0  <  x  <  X)  and  a  thawed  layer  (x  >  X).  The 
thawed  layer  is  further  divided  into  layer  2  (X  <  x  <  X  +  6)  where  temperature  changes  occur  and  layer  3  (x 
>  X  +  5)  where  thermal  effects  are  not  discernible.  We  ignore  the  finite  time  it  takes  for  the  surface  temper¬ 
ature  to  drop  to  the  freezing  point,  Tf.  This  time  will  be  small  compared  to  the  formation  time  and  a  realistic 
scenario  prior  to  the  onset  of  freezing  is  Tq  =  Tf . 


The  governing  equations  are  the  conduction  energy  equations  with  appropriate  boundary  and  initial  con¬ 
ditions;  see  the  Nomenclature  for  definition  of  symbols  not  defined  in  the  text. 

For  the  frozen  zone 


3^7;  37; 

0<x<X 

(1) 

3^2  "  dt 

7’i(X,r)  =  7’f 

(la) 

Ti  (0,  t)  =  T,. 


(lb) 


0<x<X-h5 


7 


(2) 


T2{X,t)  =  Tf 

dT2{X+h,t)_^ 

dx 

The  initial  temperature  at  the  beginning  of  freeze  i 


Tj  =  Tq  +  Gx . 


The  maximum  depth  at  any  time  to  which  the  temperature  disturbance  will  be  felt  is  X  +  8.  Then 
r2(X  +  5,r)  =  (X  +  5)G  +  7’o. 

The  energy  balance  at  the  phase  change  interface  for  the  freeze  process  is 

The  energy  balance  at  the  freezing  front  can  also  be  written  as  two  equations  (Lunardini  1981b) 


aiiM f  .  ^  mx,t)dUxA)  ^  d%{x,t) 


dx  dx 


37](X,r)  dT^jXj)  ,  ,,  \dT^{X,t) 


=  P2^«2 


d%{X,t) 


Because  of  the  initial  temperature  distribution,  during  freeze  the  heat  flow  to  the  interface  from  the 
thawed  region  will  exceed  the  geothermal  heat  flow  until  equilibrium  is  established.  Likewise,  during  a 
thaw  period  the  heat  flow  from  the  thawed  zone  will  be  less  than  the  deep  geothermal  heat  flow. 

An  approximate  solution  to  this  problem  will  be  obtained  using  the  heat  balance  integral  technique  (see 
Lunardini  1991).  In  this  method,  the  differential  equations  are  solved  on  average  over  a  finite  volume  of 
material  rather  than  at  each  point  of  the  region.  The  integration  of  the  energy  equations  over  the  regions 
where  temperature  changes  are  occurring,  0  <  x  <  X  +  5,  detailed  by  Lunardini  (1981b)  is 


X  A+0 

■^1  pici  Jt^x,  t)dx  +  P2C2  J  T2{x,t)dx  -  PifX  +  (P2C2  -  Piq)?}; 


-p2C2(x+5)  ro+Y(x+8)  |--ki 


3r|(0,r) 


Quadratic  temperature  profiles  in  regions  1  and  2  that  satisfy  the  boundary  conditions  are  chosen  as 


/  r  —  (  X  —  X 

ri=7’f  +  aiX  ^  +(«iX-Ar,)  — 
V  A  y  V  A 


72  =  7f  +  [g(5  +  2X)  +  2Ar]^^-  (gx  +  AT 


where 


a,,{AT^GX)X 

^  g  5[G(5  +  2X)  +  2Ar] 

In  general,  the  simplest  temperature  profiles  that  will  satisfy  the  boundary  conditions  should  be  chosen. 
The  accuracy  of  the  method  increases  as  the  order  of  a  polynomial  temperature  choice  increases;  however, 
the  use  of  high-order  polynomials  (third  and  higher)  is  often  not  justified  since  a  small  increase  in  accuracy 
requires  significantly  more  computational  effort.  Equation  4  can  be  used  to  find  a  relation  between  X  and  d. 
In  nondimensional  form  this  is 


Ji-fc2i[a(P  +  2)  +  2(l)]: 

g 


Equation  5,  the  energy  integral  equation,  can  now  be  written  nondimensionally,  using  eq  6  and  7  as 
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The  derivatives  of  P  and  g  can  be  found  from  the  following  equations 


<^P  p'  =  05  +<3)^4 
do  03  +  02^4 

dg  ,  0/ 


^  ^  (g  +  (|))  P(P  +  2) 

m  m 

02  =-^^^i^^[2o(P  +  l)  +  2(^] 
m 

8  *^T 


^+4  p 

g  ) 
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=  ^2l(P  +  2) 

m  =  p[a(p  +  2)  +  2(1)]. 


The  problem  has  now  been  reciuced  to  a  simple  numerical  quadrature  of  eq  9  using  the  auxiliary  relations 
of  eq  10-12.  The  numerical  solution  of  eq  9  can  be  obtained  quite  easily  with  a  personal  computer  and  any 
standard  numerical  integration  routine.  (A  FORTRAN  program,  PERM.FOR,  to  carry  out  the  integration  is 
listed  in  Appendix  E.) 

The  model  requires  only  the  ratios  of  the  thawed  to  frozen  values  of  thermal  conductivity,  specific  heat  ca¬ 
pacity  and  density  for  the  permafrost  soils.  These  property  ratios  can  be  estimated  with  good  accuracy  for  soil 
systems  as  noted  in  Appendix  A.  The  absolute  values  of  the  frozen  and  thawed  soil  properties  are  not  needed 
to  carry  out  the  solution  of  eq  9—12. 


Heterogenetic  model  verification 

It  is  possible  to  check  the  solution  for  a  special  case.  Al¬ 
though  there  is  no  exact  solution  for  the  phase-change  case,  a 
solution  was  found  for  the  transient  location  of  the  Tf  isotherm 
for  a  homogeneous  soil  with  zero  latent  heat,  i.e,  infinite  Ste¬ 
fan  number  (Lunardini,  in  prep.).  The  relation  is 


f  X  GX 

24m  Tf  -  Tg 


(13) 


Table  1.  Movement  of  Tf  isotherm,  homo¬ 
geneous  soil  for  heterogenetic  freezing. 

St  =1000, £  =  0,Tf  -T=  10°C, 

(j)  =  0,  G  =  0.0286“C/m. 


X 

(m) 

Time  (years) 

Ed  9  Exact  eq  13 

%  difference 

27.97 

2.08 

2.15 

3.6 

69.93 

22.63 

25.30 

10.6 

139.86 

204.48 

234.73 

12.9 

314.69 

53,053.2 

48,141.3 

-10.2 

332.17 

242,852.4 

238,469.1 

-1.8 

If  we  let  the  Stefan  number  be  large  and  hold  the  property 
ratios  to  unity,  the  approximate  solution  can  be  compared  to  this  exact  relation.  Table  1  notes  the  results  for  a 


typical  case. 

The  comparison  indicates  that  the  approximate  technique  gives  good  results,  especially  as  the  time  in¬ 
creases.  The  results  also  show  that  the  Tf  isotherm  requires  surprisingly  long  times  to  penetrate  deeply,  even 
without  phase  change.  This  is  explainable  by  the  sensible-to-latent  heat  ratios  to  be  examined  below. 


Sensible  and  latent  heats 

The  total  energy  extracted  from  a  unit  area  of  soil  is  the  sum  of  the  latent  and  sensible  energies. 
<2t  =  2l  +  Gs  • 


(14) 


The  latent  energy  is 

eL  =  XL  (15) 

while  the  sensible  heat  flow  is 

X  X  x+b 

Qs  =  Cfj  [Tf  -  7’,(x,/)Vx  +  cj(7]  -  Tf)  dx  +  C,  j 

()  0  ^ 

where  Ti  (x)  is  the  original  temperature  before  the  freeze 
leads  to 

Qs  =  CfAT,  XIC2J  [(a/2)  +  ^  +  P(a-i-  ([))  /  3]  -i-  (0.5/g  +  1)  /  3}  .  (17) 

The  ratio  of  the  sensible  to  the  latent  heat  is 

Qs  IQh  =  [Ql  [(f^(2)  +  ([)  +  P(a  +  ([))/  3]  +  (0.5/g  +  l)/3}.  (18) 

This  ratio  is  quite  large  even  for  small  Stefan  numbers  and  tends  to  increase  as  the  freeze  depth  increases. 


[T^-T4x,t)]dx  (16) 

starts.  Using  the  temperature  relations  (eq  6  and  7) 
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Neumann  solution 

The  Neumann  solution  is  a  special  case  that  occurs  when  the  initial  temperature  gradient  G  is  zero. 
Thus,  it  is  the  conduction  solution  that  will  predict  the  minimum  time  needed  to  form  a  given  thickness  of 
frozen  material.  The  exact  solution  is  well  known  as 

X  =  2YnV^. 

Yn  is  a  function  of  ,  (j)  and  the  property  ratios  and  is  known  (Lunardini  1981).  Equation  19  can  also  be 
written  as 


or 


(23) 

(24) 


Effect  of  quaternary  freeze-thaw  cycles 

The  previous  discussion  assumed  that  the  soil  was  initially  unfrozen,  with  a  temperature  distribution  giv¬ 
en  by  eq  2c,  as  shown  in  Figure  9.  However,  the  cyclic  warming  and  cooling  that  has  taken  place  over  the 
past  million  years  will  tend  to  cyclically  lower  the  ground  temperature  as  compared  to  the  previous  initial 
temperature  distribution.  If  the  temperature  is  lower  than  that  of  eq  2c,  the  permafrost  growth  will  be  accel¬ 
erated.  Figure  10  shows  the  thermal  state  of  permafrost  at  equilibrium  with  a  geothermal  heat  flow — the 
solid  line.  The  approach  to  the  equilibrium  state  of  the  per¬ 
mafrost  is  affected  by  the  initial  temperature  of  the  soil  in 
that  the  time  to  reach  equilibrium  will  be  greatly  changed 
but  the  final  equilibrium  thickness  is  not  affected  by  the 
initial  assumption  (see  eq  23). 

Let  us  assume  that  surface  warming  to  7f  occurs  and 
the  permafrost  eventually  reaches  a  stage  where  the  entire 
frozen  zone  is  at  7f .  The  system  starts  to  thaw  in  this  con¬ 
dition  and  we  ignore  the  time  needed  to  reach  this  state 
(small  compared  to  melt  time  scales)  and  also  neglect  the 
bottom  freeze-thaw  during  this  initial  warming  time.  Fig¬ 
ure  1 1  shows  the  thermal  state  as  the  bottom  thaw  com¬ 
mences,  at  X  =  Xq.  The  geothermal  heat  flow  is  used  to  Figure  10.  Temperature  distribution  after 
melt  permafrost  and  also  to  warm  the  soil  in  the  region  X  initial  freeze  of  soil. 


long-term  equilibrium  (b). 


+  8.  The  maximum  thaw  rate  occurs  if  all  of  the  geothermal  heat  flow  goes  into  melting,  but  this  is  physi¬ 
cally  impossible. 

Given  sufficient  time,  the  entire  permafrost  volume  melts  and  the  soil  temperature  is  as  shown  in  Figure 
12— curve  a.  A  layer  of  soil  of  thickness  +  8o  will  be  thermally  modulated.  The  temperature  gradient  for 
depths  greater  than  Xo+  5o  is  G,  the  geothermal  gradient.  Also  shown  in  Figure  12  is  the  dashed  line  (curve 
b)  denoting  the  equilibrium  temperature  distribution  with  the  same  geothermal  gradient  but  no  cooling. 
This  state  would  be  reached  if  the  surface  remained  at  Tf  for  a  long  time  after  thaw  was  completed.  Clearly, 
the  sensible  heat  that  must  be  removed  in  the  modulated  case  (curve  a)  will  be  less  than  that  for  the  original 
freeze  discussed  earlier  (curve  b).  This  means  that  cooled  soil  (after  cyclic  melting  or  warming)  will  freeze 
much  faster  with  the  same  mean  surface  temperature.  This  temperature  modulation  will  only  be  significant 
for  zones  of  relatively  shallow  permafrost,  for  which  appreciable  thaw  can  take  place  during  the  intergla¬ 
cials  of  about  15,000  years.  The  equations  and  solution  are  discussed  in  Appendix  B. 


Syngenetic  growth  of  permafrost 

Syngenetic  growth  of  permafrost  oc¬ 
curs  when  material  is  deposited  at  the 
surface  while  freezing  is  in  progress. 
This  is  inherently  a  much  more  efficient 
freeze  process  and  the  growth  of  frozen 
layers  can  be  greatly  accelerated.  Figure 
13  shows  a  sketch  of  the  process  with  the 
surface  deposit  laid  down  such  that  the 
surface,  held  at  a  constant  temperature 
Ts,  is  moving  at  a  constant  velocity  U. 
The  total  frozen  zone  at  any  time  is  equal 
to  the  material  frozen  at  the  interface 
X(0,  plus  the  depositional  layer  U^.  The 
syngenetic  system  will  be  inherently  un¬ 
stable  with  no  equilibrium  solution.  If 
the  frozen  zone  should  equal  the  hetero- 
genetic  steady-state  value  (JJ  =  0),  the 
motion  of  the  upper  surface  will  cause 
melting  at  the  base  of  the  permafrost. 


Figure  13.  Geometry  and  coordinate  system  for  freeze  of  a 
semi-infinite  medium  with  moving  upper  surface. 
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Thus,  at  some  time  during  the  formation  process,  the  bottom  must  be  melting  to  compensate  for  the  deposi- 
tional  growth  but  the  process  cannot  be  steady. 

Let  us  fix  the  upper  surface  so  that  it  remains  stationary,  as  shown  in  Figure  14.  Then  it  would  appear 
that  a  steady  flow  of  material  is  moving  at  a  constant  velocity  U  and  the  original  soil  surface  seems  to  be 
moving  downward  at  a  steady  velocity.  The  energy  equation  is  transformed  such  that  a  convective  term  ex¬ 
ists.  The  equations  for  regions  1  and  2  are 


a 


dx^ 


■  u 


dTj  dTj 

dx  dt 


=  0 


^’  =  1,2. 


(25) 


The  boundary  and  initial  conditions  are  exactly  the  same  as  those  of  the  heterogenetic  case  (see  Appendix 
C  for  details  of  the  syngenetic  equations). 

The  heat  balance  integral  form  of  the  energy  equations  is 


— {Pl^i  j  71  {x, t)dx  +  P2C2  j  T2  {x,  t)dx  -  pi«  +  (P2C2  -  PiC]  )rf X 

0  X 


-P2C2(X  +  8) 


=  -h 


+  k2G  -  Picif/Ari  -  P2C2t/f AT  +  G(X  +  8)1. 

dx 


(26) 


We  note  that  this  is  identical  to  the  relation  for  heterogenetic  growth,  eq  5,  except  for  the  two  additional 
terms  on  the  right-hand  side  of  eq  26.  Carrying  out  the  integrations  and  making  eq  26  nondimensional, 
leads  to  the  following  result. 


t  = 


a 


jKida 


0 


+^>2^-  — 
og 

f.-sl) 

g  J 

-C2\0\ 

1  1 
a’ 

fi-2l 

j 

+  k2i-\V 

|l  + 1  /  ♦S'y  +  C21 1 

[(^  +  a(p  +  l)]} 

where 


(27) 


(28) 
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\|/  = - ^  . 

^  Gaj 

Now  the  integration  of  eq  27  follows  exactly  as  was  done  previously  (see  FORTRAN  program 
PFTS YNB .FOR  in  Appendix  E).  The  model  requires  only  the  ratios  of  the  thawed  to  frozen  values  of  ther¬ 
mal  conductivity,  specific  heat  capacity,  and  density  for  the  permafrost  soils,  as  was  noted  earlier. 

Syngenetic  model  verification 

It  is  possible  to  check  the  solution  for  a  special  case  as  was  done  for  the  heterogenetic  growth.  Although 
there  is  no  exact  solution  for  the  phase-change  case,  there  is  an  exact  solution  for  the  transient  location  of 
the  Tf  isotherm  for  the  same  problem  with  a  homogeneous  soil  with  zero  latent  heat,  i.e,  infinite  Stefan 
number  (Lunardini,  in  prep.).  The  relation  is 

gWf(5_i)  ei-fc— ^-(A-hl)  erfc— ^  +  2A  =  0  (29) 

2'V'c  l-jx 


where  A  =  Gf  ~  \j/x 
5  =  Gf  +  Vj/T 

Gf  =  location  of  7f  isotherm. 

If  we  let  the  Stefan  number  be  large  and  hold  the  property  ratios  to  unity,  the  heat  balance  integral  solution 
for  syngenetic  growth  can  be  compared  to  this  exact  relation.  Table  2  notes  the  results  for  typical  cases. 


Table  2.  Movement  of  Tf  isotherm,  homo¬ 
geneous  soil,  syngenetic  freezing. 

St  =1000,  ^  =  0,  Tf-  Ts  =  10°C,  (j)  =  0,  G  = 
0.0286'^C/m. 


T 

Gf  exact 

Of  approximate 

%  difference 

a. 

U=\0  mm/yr 

0.0014 

0.0087 

0.090 

-1.5 

0.0069 

0.1634 

0.1700 

^.0 

2.033 

0.7739 

0.7700 

0.5 

2.8232 

0.8206 

0.8100 

1.3 

b. 

U=1  mm/yr 

0.005 

0.1448 

0.15 

-3.6 

0.099 

0.3874 

0.4006 

-3.42 

1.0689 

0.6572 

0.6662 

-1.37 

11.4705 

0.8795 

0.8720 

0.85 

The  results  indicate  that  the  approximate  technique  gives  excellent  results,  especially  as  the  time  in¬ 
creases.  Thus,  the  Heat  Balance  Integral  method  and  the  numerical  quadrature  are  robust  even  for  very  long 
time  spans.  See  also  eq  C17  for  further  verification  of  the  solution  method  with  phase  change. 

DISCUSSION 

Equation  9  was  solved  numerically  using  Simpson’s  rule.  This  resulted  in  values  of  the  permafrost  depth 
versus  time  as  a  function  of  Sj,  e  and  the  thermal  property  ratios  of  the  frozen  and  the  thawed  zones.  The 
results  are  presented  in  Figures  15-17. 

These  graphs  depend  only  upon  the  quantities  ^t,  <|)  and  £.  In  Appendix  A  it  is  shown  that  the  soil  por¬ 
osity,  e,  determines  the  saturated  soil  property  ratios.  The  thermal  property  ratios  used  for  the  graphs  are 
listed  in  Tables  A1  and  A2.  The  graphs  are  only  valid  for  the  particular  soil  ratios  given.  However,  this 
does  not  affect  the  validity  of  the  model  itself  Any  specific  site  can  be, modeled  by  using  site-specific  prop¬ 
erty  ratios  in  eq  9.  Figures  15-17  can  be  used  to  estimate  permafrost  formation  times  for  a  wide  range  of 
surface  temperatures  and  geothermal  gradients. 

The  graphs  can  also  be  applied  to  variable  surface  temperatures  with  a  bit  of  manipulation.  Figure  18 
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Figure  15.  Formation  time  of  permafrost,  saturated  mineral  soils,  (|)  =  0,  e  =  0.4. 


Figure  16.  Formation  time  of  permafrost,  saturated  mineral  soils,  (])  =  0,  E  =  0.379  (Prudhoe  Bay,  Alaska). 

plots  the  time  needed  to  reach  90%  of  the  equilibrium  permafrost  thickness  at  a  site.  Figure  19  shows  the 
thickness  of  permafrost  formed  after  15,000  years  as  a  function  of  average  surface  temperature  for  different 
soils  as  characterized  by  the  soil  porosity.  Figures  20  and  21  show  that  syngenetic  growth  greatly  reduces 
the  formation  time  during  the  latter  stages  of  permafrost  growth.  However,  there  may  be  significant  ques¬ 
tions  as  to  whether  surface  deposition  has  continued  for  long  periods  of  time  at  any  given  location. 
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Figure  17.  Formation  time  of  permafrost,  saturated  mineral  soils,  (j)  =  0,  8  =  (9.5. 


Figure  18.  Formation  time  for  permafrost  to 
reach  90%  equilibrium  thickness,  saturated 
mineral  soils. 


Figure  19.  Effect  of  soil  porosity  on  permafrost 
thickness  after  15,000  years,  saturated  mineral 
soils,  G  =  0.0286°C/m. 


Figure  20.  Syngenetic  growth  of  permafrost,  saturated  mineral 
soils,  (j)  =  0,  e  =  0.379,  Sp  =  0.15  (Prudhoe  Bay,  Alaska). 


Figure  21.  Syngenetic  growth  of  permafrost,  saturated  mineral 
soils,  (j)  =  0,  £  =  0.379,  \}/  =  0.01  (Prudhoe  Bay,  Alaska). 


Table  3.  Paleotemperature  scenarios,  Prudhoe  Bay,  Alaska 
(after  Osterkamp  and  Gosink  1991). 


Scenario 

T  * 

*.v 

(^C) 

Duration  before 
present  time 
(years) 

Stefan 

number 

Present  (Lachenbruch  et  al.  1982) 

-10.99 

10,000-15,000 

0.1438 

Brigham  and  Miller  (1983) 

-13.69 

240,000 

0.1827 

Matteucci  (1989) 

-11 

300,000 

0.1440 

Vostok,  Robin  (1983)+ 

-16 

160,000 

0.2159 

E.  Siberia  (Maximova  and 
Romanovsky  1988)+ 

-11.3 

170,000 

0.1483 

*  Average  value  of  the  fluctuating  surface  temperature  over  the  indicated 
time  span. 

t  Fitted  to  present  Prudhoe  Bay  surface  temperature. 
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Prudhoe  Bay,  Alaska 

Considerable  information  on  the  permafrost  is  available  from  oil  wells  in  the  Prudhoe  Bay,  Alaska,  area 
(Lachenbruch  et  al.  1982).  Using  the  actual  permafrost  data,  we  note  that  the  property  ratios  for  e  =  0.379 
(Tables  Al  and  A2)  are  very  close  to  measured  and  estimated  values  (Lachenbruch  et  al.  1982).  Some  tem¬ 
perature  possibilities  for  Prudhoe  are  listed  in  Table  3.  Figures  16a  and  b  can  be  used  to  estimate  the  per¬ 
mafrost  formation  time,  depending  upon  the  temperature  chosen. 

Example  I 

At  Prudhoe  Bay,  Alaska,  the  permafrost  has  the  following  present  conditions: 


-10.99°C 

G  = 

0.0286°C/m 

e  = 

0.379 

Sj  = 

0.1440 

ai  = 

58.89  m^/yr 

kg  - 

4.34  W/m  K 

measured  k2\  = 

0.5795 

measured  present  permafrost  thickness,  Xp  = 

599.3  m. 

Then,  the  equilibrium  permafrost  thickness  is 

X.  =-^  =  602.8  m. 

hxG 


The  calculated  equilibrium  thickness  is  essentially  the  same  as  the  measured  value.  How  long  would  it  take 
to  reach  this  depth  if  the  surface  temperature  had  been  =  -10.99°C  for  an  indefinite  period?  The  nature 
of  the  solution  is  such  that  the  final  equilibrium  values  will  be  reached  only  after  very  long  times.  From 
Figure  1 8  we  find  the  time  to  reach  90%  of  equilibrium  (542.0  m)  is 


:ai 


r  =  241.2 


thus  t  =  500,740  years.  This  time  is  obviously  quite  long  and  suggests  that  the  present  climate  of  Prudhoe 
Bay  is  probably  considerably  warmer  than  it  has  been  on  average  over  the  past  glacial  cycles.  Such  warm¬ 
ing  over  the  past  15,000  years  is  widely  accepted. 


Example  2 

Prudhoe  Bay,  with  the  Brigham  and  Miller  (1983)  paleotemperature  scenario,  has  the  following:  Tj  =  - 
13.69°C  for  225,000  years  before  warming  to  -1 1°C  in  the  last  15,000  years.  For  this  case  Xg  =  763.5  m. 
For  t  =  225,000  years,  X  =  67.3.  From  Figure  16,  with  Si  =  0.1827,  we  read  a  =  1.412.  Thus,  in  225,000 
years  the  permafrost  will  grow  to  X  =  626.5  m.  This  value  will  then  slowly  decay  to  the  new  equilibrium  of 
602.8  m  over  15,000  years.  This  requires  a  melt  rate  of  1.58  mm/yr.  This  value  is  well  within  the  estimated 
average  thaw  rate  of  2.5  mm/year  for  this  case  (see  Appendix  D).  Note  that  a  lower  surface  temperature 
will  greatly  accelerate  growth  since  the  new  equilibrium  depth  will  be  greater  than  before.  Hence,  a  much 
larger  fraction  of  the  growth  will  be  during  the  early,  rapid  growth  stage. 

Figure  22  shows  the  permafrost  thickness  at  Prudhoe  Bay  after  six  glacial  cycles  with  some  typical  fea¬ 
tures  of  permafrost  growth  demonstrated.  First,  the  initial  permafrost  growth  is  quite  rapid,  reaching  a 
thickness  of  570  m  after  120,000  years,  with  the  paleotemperature  model  of  Figure  7a.  The  thickness  then 
slowly  approaches  an  equilibrium  value  of  739  m  but  it  will  surpass  the  present  thickness  of  600  m  after 
about  185,000  years.  Thus,  a  paleotemperature  model  as  cold  as  that  of  Figure  7b  will  yield  permafrost  that 
is  too  thick.  Also  shown  is  the  finite  difference  prediction  of  Osterkamp  and  Gosink  (1991),  using  Figure  6 
(which  has  the  same  mean  temperature  as  Fig.  7a)  and  starting  with  600  m  of  permafrost,  that  indicates 
much  faster  permafrost  growth  and  thicker  permafrost.  Their  quasi-steady  model  neglects  sensible  heat,  as- 
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Figure  22.  Growth  of  permafrost  at  Prudoe  Bay,  Alaska;  effect  of  paleotemperature  scenarios. 


sumes  fixed  geothermal  heat  flow  to  the  freezing  inter- 
face,  and  starts  with  permafrost  depth  far  less  than  equi¬ 
librium.  All  of  these  approximations  overestimate  the 
freeze  rate.  The  second  curve  of  Figure  22  shows  the  sig¬ 
nificant  effects  if  the  paleotemperature  model  is  modified 
by  only  a  small  amount.  The  predicted  permafrost  thick¬ 
ness  will  reach  present  values  after  about  640,000  years 
and  will  tend  to  oscillate  about  this  value. 


Table  4.  Effect  of  previous  cooling  on  perma 
frost  growth  time.*  X  =  541.0  m. 


Case 

G 

Time 

(years) 

Comment 

No  pre-cool 

0.0286 

486,100 

uses  eq  10 

Pre-cool 

0.0220 

78,900 

uses  eq  10 

Pre-cool 

0.0286 

105,400 

uses  eq  B15 

Syngenetic 

0.0286 

172,500 

eq  21,U  mm/yr 

*  Soil  properties  for  Prudhoe  Bay,  ATi  -  10.0°C. 


Example  3 

Consider  the  effect  of  previous  cooling  at  Table  5.  Formation  time  of  deep  permafrost.* 


Prudhoe  Bay.  If  the  equivalent  geothermal  gra¬ 
dient  is  G  =  0.0220°C,  then  the  time  to  reach 
541  m  is  78,940  years.  The  results  of  three  pos¬ 
sibilities  are  shown  in  Table  4.  Note  the  very 
large  effect  of  previous  cooling  or  cyclic  ther¬ 
mal  modulation. 

Deep  permafrost 


=  1813  m. 


Time 


Case 

G 

(years) 

Comment 

No  pre-cool 

0.0286 

4,190,600 

uses  eq  10 

Pre-cool 

0.0220 

488,900 

uses  eq  10 

Pre-cool 

0.0286 

358,900 

uses  eq  B 1 5 

Syngenetic 

0.0286 

784,100 

eq  27,  U  =  1  mm/yr 

Neumann  ((])  =  0) 

0.0 

64,800 

absolute  minimum  time 

Example  1  *  properties  for  Prudhoe  Bay,  AT\  -  29.3°C. 

Consider  the  case  of  very  thick  permafrost, 

on  the  order  of  1600  m.  Let  the  properties  be  those  of  Prudhoe  Bay  but  e  =  0.4,  AT\  =  29.27°C,  G  =  0.0286, 


ai  =  57.99  m^/yr,  then  Xq  =1813  m.  The  value  of  the  surface  temperature  chosen  is  on  the  order  of  12°C 


less  than  present  winter  temperatures  experienced  in  parts  of  Canada,  Russia  and  Greenland,  although  it  is 
doubtful  that  such  temperatures  could  have  persisted  for  1  million  years.  The  value  used  illustrates  the  long 
time  needed  to  form  deep  permafrost  by  conduction  alone.  We  will  find  the  time  required  to  form  90%  of 
the  permafrost  or  1632  m.  The  calculations  are  as  before,  with  Sj  =  0.4  (see  Fig.  15  for  the  case  without 
previous  cooling).  The  results  are  given  in  Table  5.  The  formation  time  is  very  long,  even  with  previous 
cooling.  This  example  used  Prudhoe  Bay  properties  and  geothermal  gradient,  which  could  be  significantly 
different  at  a  site  in  Siberia  with  deep  permafrost. 


Example  2 

Consider  the  question  of  the  maximum  permafrost  thickness  that  is  probable.  This  will  occur  if  the  fro¬ 
zen  thermal  conductivity  is  large,  the  geothermal  heat  flow  and  latent  heat  are  low,  and  the  surface  temper¬ 
ature  is  minimal  (within  the  constraints  discussed  earlier).  Let  e  =  0.2,  =  -23.5°C,  =  0.042  W/m^,  and 

assume  coarse-grained  soil  with  =  5.86  W/m  K.  The  results  of  the  calculation  are  shown  in  Table  6. 
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The  permafrost  reaches  a  thickness  of  2132  m  after  1  million 
years  and  a  value  of  2255  m  after  2  million  years.  Balobaev  et  al. 
(1978)  note  that  the  greatest  permafrost  thicknesses  recorded  are  on 
the  East  Siberian  platform  and  present  graphs  with  maximum  perma¬ 
frost  thicknesses  of  1500  m.  These  thicknesses  are  on  the  order  of  the 
value  predicted  here  and  it  is  not  likely  permafrost  much  thicker  than 
this  has  ever  existed,  since  the  required  time  exceeds  the  plausible 
time  available.  These  extreme  thicknesses  are  not  in  thermal  equilib¬ 
rium  with  the  present  surface  temperatures  and  are  slowly  thawing. 


Table  6.  Extreme  predicted  perma¬ 
frost  thickness. 


X(m) 


Time 

(years) 

X(m) 

hetero  genetic 

syngenetic 

U  -  /  mm/yr 

100,000 

1488.4 

1539.9 

532,257 

2000.0 

2158.6 

1,000,000 

2132.0 

2362.9 

2,000,000 

2254.5 

2529.5 

Infinity 

2600.2 

2600.2 

CONCLUSIONS 

The  calculations  and  examples  indicate  that  the  growth  of  permafrost,  with  pure  conduction  heat  trans¬ 
fer,  is  governed  by  the  transient  surface  temperature,  the  geothermal  heat  flow  and  the  soil  thermal  proper¬ 
ties.  Permafrost  grows  very  rapidly  for  an  initial  phase  and  then  asymptotically  approaches  a  steady-state 
value  after  time  spans  of  immense  length.  Very  thick  permafrost  may  have  required  the  total  Quaternary 
Period  to  form.  It  is  likely  that  permafrost  is  not  in  equilibrium  at  most  sites.  The  bottom  growth  and  decay 
of  permafrost  are  so  slow  that  accurate  methods  of  detecting  which  is  occurring  (or  if  equilibrium  exists) 
are  not  available  for  the  field.  Permafrost  less  than  600  m  can  grow  within  50,000  years,  with  surface  tem¬ 
peratures  only  slightly  lower  than  present  values,  but  deeper  permafrost  depths  require  time  scales  of  sever¬ 
al  ice  ages  and  quite  low  temperatures  to  form. 
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APPENDIX  A:  SOIL  PROPERTIES  AND  RATIOS 


The  thermal  conductivity  of  a  mixture  such  as  a  soil  can  be  estimated  by  using  the  weighted  geometric 
mean  (Lachenbruch  et  al.  1982,  Gold  and  Lachenbruch  1973,  Lunardini  1981a).  This  can  be  written  for  a 
general  soil  as 

k={k^y^{k^Y'“{KY^{kiY''  (Al) 

where  kg,  ki  and  /:w  are  the  thermal  conductivities  of  air,  soil  solids,  ice  and  water;  x^,  Xg,  and  are 
the  volumetric  fractions  of  air,  soil  solids,  ice  and  water.  The  geometric  mean  is  usually  better  than  the  as¬ 
sumption  of  parallel  geometry  (weighted  arithmetic  mean),  which  is  often  used  for  simplicity. 

Saturated  soil 

Many  assumptions  can  be  made  concerning  the  soil  saturation  and  porosity  but  simple  approximations 
will  be  used  here.  If  the  soil  is  always  saturated,  has  a  constant  void  ratio  £,  and  all  of  the  water  freezes, 
then  Xyv  =  e  and  the  conductivity  ratio  can  be  expressed  as 


(A2) 

(A3) 

(A4) 

where  y  =  0.9825  is  a  temperature  correction  for  kg  (Lachenbruch  et  al.  1982)  and  =  1.34/5.45  = 
0.2459. 

The  volumetric  specific  heat  for  the  system  may  be  expressed  as  follows,  for  the  thawed  and  frozen 
states 


“■  ^SU  (1  "”^w)  ^w-^w 

Cf  =  Csf(l-Xi)  +CiXi  (A6) 


where  Cgu  and  Csf  are  the  volumetric  specific  heats  of  unfrozen  and  frozen  solids,  and  and  Q  are  the 
volumetric  specific  heats  of  water  and  ice. 

It  is  fortunate  that  the  volumetric  specific  heats  of  soil  solids  and  ice  are  all  about  the  same.  For  exam¬ 
ple,  the  specific  heat  of  organic  solids  is  0.461  cal/cm^  °C,  for  mineral  solids  it  is  0.420,  and  for  ice  it  is 
0.459  (Lunardini  1981a).  If  one  assumes  that  the  values  for  the  solids,  except  for  ice,  change  little  through 
the  phase  change  then 


Cf  =  0.4202  +  0.03888 

C^_C  0.4296 +  0.5708£ 
~Q~  ”  0.4202  +  0.0388e ' 

The  density  ratio  is 

PiL-  l-0.6154e 

Pf  ~  l-0.650e  ■ 


(A7) 

(A8) 


(A9) 
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Finally,  the  ratio  of  thermal  diffusivities  is 


a„ 


“f 


=  «21 


■  «21 
C21 


The  latent  heat  is 


L  =  79.71  8 
The  Stefan  Number  is 


(AlO) 


(All) 


„  C,,  ,  C,A71  .  0.4202  +  0.0388S 

Sr-yW-r.) - -  ATI.  (A12) 

It  is  possible  to  present  the  results  for  soil  systems,  quite  efficiently,  since  the  property  ratios  can  be  de¬ 
scribed  as  functions  of  the  soil  void  ratio  8  (Lunardini  and  Varotta  1981).  Using  the  thermal  conductivities 
of  Table  Al,  the  property  ratios  used  in  the  calculations  are  given  in  Table  A2.  The  thermal  conductivity 
ratio  will  be  representative  of  soil  that  is  not  too  dry.  Thus,  eq  A4  and  AlO  should  be  acceptable  if  8  >  0.2 
(Kersten  1949). 


Table  Al,  Thermal  conductivity  of 
materials. 

Thermal  conductivity 


Substance _ W/(m°C) 


Water  0.561 

ice  2.281 

air  0.0237 

silicaceous  soil  solids*  4.29-5.87 


*Lachenbruch  et  al.  (1982) 

Balobaev  et  al.  (1978)  note  that  for  limestone 
-  0.021  W/m^  at  60-800  m,  anomalously  low  heat 


Table  A2.  Calculated  saturated  granular 
soil  thaw-freeze  property  ratios. 


Soil 

porosity  e 

eq  21 

CJCf 
eq  24 

ajuf 
eq  26 

pjpf 
eq  25 

0.2 

0.7448 

1.2706 

0.5862 

1.008 

0.3 

0.6484 

1.3909 

0.4662 

1.0129 

0.379 

0.5812 

1.4847 

0.3915 

1.0174 

0.4 

0.5645 

1.5094 

0.3740 

1.0187 

0.5 

0.4915 

1.6265 

0.3022 

1.0256 

dolomite,  k  =  2.44  -  3.37  W/m  °C  and  =  0.017 
values. 


Nonsaturated  soil 

For  the  nonsaturated  soil,  assuming  that  the  porosity  does  not  change  during  phase  changes,  the  ratio  of 
thawed  to  frozen  thermal  conductivity 


(A13) 

(A14) 

(A15) 

where  S  is  the  thawed  soil  saturation  level,  p^i  =  Pw  /pi;  Pi^Pw  “  0-91.  Interestingly,  the  ratio  k2i  can  have 
the  same  values  as  for  the  saturated  case  if  the  saturation  has  certain  values,  e.g.,  S  =  0.756,  8  =  0.379. 
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APPENDIX  B:  QUATERNARY  CYCLIC  THERMAL  MODULATION 


Thaw  process 

Assuming  that  melting  occurs  with  a  fraction /of  the  geothermal  heat  flow  going  into  increasing  the  sen¬ 
sible  heat  allows  an  estimate  to  be  made  of  the  thaw  time  and  the  final  temperature,  noted  in  Figure  12, 
which  will  then  be  the  initial  temperature  distribution  for  the  next  freeze  cycle.  Referring  to  Figure  1 1,  we 
assume  that  the  fraction  of  geothermal  energy  available  for  melt  at  any  time  is 


(/^-l)X-HXo 


(Bl) 


where /n,  is  the  melt  fraction  at  the  conclusion  of  melt  and  Xq  is  the  initial  frozen  thickness.  The  tempera¬ 
ture  in  the  region  of  changing  temperature  is 

T2  =  T{+  b{x  -  Ai)  +  c{x  -Ai)2  (B2) 


where  Ai  =  Xq  -X  and  b,c  are  functions  of  time  only.  The  initial  temperature  of  the  soil  is 


fTf  x<Xo 

[T(+G{x-X^)  x>X^' 


(B3) 


The  sensible  heat  added  to  the  soil  is 


X„  Xo+5 

=  j  {T^-TCjdx  +  C^  I  {T2-Ti)dx. 

Xo-X  Xo 

The  sensible  heat  addition  is  shown  as  the  shaded  region  of  Figure 
Bl.  The  sensible  heat  at  any  instant  is  then 

a=iGC„x(5-|].  (B5) 

The  change  in  sensible  heat  must  equal  the  fraction  of  the  geother¬ 
mal  energy  not  used  for  melting,  or 

^  =  (B6) 

The  energy  balance  at  the  melt  interface  is 

L—  =  fq„.  (B7) 

dt  ® 

The  solutions  to  these  equations  lead  to 


(B4) 


Figure  Bl.  Bottom  melt  process. 


1  f 
-GQX  5 
3  “  I 


(B8) 


X  = 


~  fm) 


(B9) 
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where 


^^(/m  l)gg 

X^L 

The  time  to  complete  melt  is 


_  XpLln/^ 

This  leads  to  the  value  of  5  when  melt  is  completed. 


(BIO) 


5. 


2  GC„ 


v/m-1 


•1 


(BID 


We  define  a  linear  temperature  distribution  that  will  have  the  same  sensible  heat  when  thaw  is  completed. 


7;'=7f +G'a: 


G  5„+X/ 


(B12) 

(B13) 


where  G’  is  the  equivalent  geothermal  gradient.  Finally,  the  new  temperature  distribution  at  the  beginning 
of  freeze  is 


7j  =  Tf  +bf^x  +  CQX^ 


(B14) 


where 

G(5,-X,)  GX, 

5o  +  X,  '  °  (5„+X„f 

This  initial  temperature  distribution  is  shown  as  curve  a 
in  Figure  12.  Table  B1  shows  some  results  for  Prudhoe 
Bay.  Note  the  long  melt  times  even  if /is  as  high  as 
90%. 

Freeze  of  cooled  soil 

The  freezing  process  is  as  discussed  earlier  except 
that  the  initial  soil  temperature  is  lowered  as  noted  in 


Table  Bl.  Melt  relations  for  Prudhoe  Bay. 


f 

f,  „  , 

8o 

(in) 

G' 

G 

Melt  time 
(years) 

f ,  -1 

0.1 

0.55 

1.5584 

7948.1 

0,9064 

108,926 

0.2 

0.60 

1.0118 

5265.6 

0,8636 

85,654 

0.3 

0.65 

0.720 

3833.5 

0.8196 

73,230 

0.4 

0.70 

0.5272 

2887.3 

0.7706 

65,022 

0.5 

0.75 

0.3863 

2195.8 

0.7139 

59,023 

0.6 

0.80 

0.2771 

1659.9 

0.646 

54,373 

0.8 

0.90 

0.1157 

867.8 

0.455 

47,502 

Xo  =  600  m,  L  =  30.21  cal/cm^,  =  0.6457  cal/cm^  °C. 
=  1.35  X  10-6  cal/cm2  s,  G  =  0.0286°C. 


Figure  10.  The  basic  equations  used  earlier  are  still  valid 

except  that  the  coefficients  of  eq  6  and  7  change,  owing  to  the  new  initial  temperature  given  by  eq  B14. 


The  basic  equation,  replacing  eq  9,  is 


dFa 

dx 


0 

8J 


(B15) 


1  I  P2I 
3  Sj  j 


-C2,0 


3  2 


-.^(P-i)(P+iy 


(B16) 


P 

8 


^2ia[-p  +  2mo(p  +  l)  +  2a5o(p  +  lf 


2p2l(g-l)P 

5t 


(B17) 
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(B18) 


«2i(6-P) 

P(2e-P) 

where  Q  -  mo  (P+1)  +  oSq  (p+1)^ 
mo  =  bo/G 
Sq  =  cqATj/G^. 

The  solution  to  these  equations  follows  in  exactly  the  same  way  as  previously  discussed. 
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APPENDIX  C:  HEAT  BALANCE  INTEGRAL  EQUATIONS 
FOR  SYNGENETIC  PERMAFROST  GROWTH 


Method  1 

We  will  formulate  the  basic  equations  in  terms  of  a  convective  system  with  mass  flowing  through  the 
stationary  upper  surface  at  constant  velocity  U,  as  shown  in  Figure  14.  The  governing  equations  are 

a.^-U^-^  =  0  0<x<X  (Cl) 

‘  dx^  dx  dt 


7’i(X,r)  =  7f  (Cla) 

7i(0,0  =  Ts  (Clb) 


«2 


dx^  dx  dt 


X<x<X+d 


(C2) 


T2(X,t)  =  Tf  (C2a) 

dT2iX  +  5,t) 
dx 


r2(X  +  8,r)  =  (X  +  5)G  +  7’o.  (C2c) 

The  initial  temperature  at  the  beginning  of  freeze  is 

Ti=Tg+  Gx.  (C2d) 

The  energy  balance  at  the  phase  change  interface  for  the  freeze  process  is 


dx 


iX,t)-b 


dT2 


dX 


dx  dt 


(C3) 


Because  of  the  initial  temperature  distribution,  during  freeze,  the  heat  flow  to  the  interface  from  the  thawed 
region  will  exceed  the  geothermal  heat  flow  until  equilibrium  is  established.  Likewise,  during  a  thaw  per¬ 
iod,  the  heat  flow  from  the  thawed  zone  will  be  less  than  the  deep  geothermal  heat  flow.  The  energy  bal¬ 
ance  at  the  freezing  front  can  also  be  written  as  two  equations  (Lunardini  1981b) 


-k 

-k 


dT^jXj) 

dx 


+  k 


dT2iX,t)  dTi(X,t) 


dx 


dx 


dT,{X,t)  dT^jXj)  ^  , 
dx  dx 


dT^jXj) 

dx 


=  Pl^«l 


d^T]{X,t) 

dx^ 


=  P2'(«2 


d%(X,t) 

dx^ 


(C4) 

(C5) 


Quadratic  temperature  profiles  in  regions  1  and  2  that  satisfy  the  boundary  conditions  are 

r,  =rf +  (C6) 


29 


T2=Tf+  [G(6  +  2X)  +  2AT]  f^-^1  -  {GX  +  AT) 


{Cl) 


where 


ayX- 


AT 


a^,{AT  +  GX)X 
^  5[G(8  +  2X)  +  2Ar] 


In  general,  the  simplest  temperature  profiles  that  will  satisfy  the  boundary  conditions  should  be  chosen. 
The  accuracy  of  the  method  increases  as  the  order  of  a  polynomial  temperature  choice  increases;  however, 
the  use  of  high-order  polynomials  (third  and  higher)  is  often  not  justified  since  a  small  increase  in  accuracy 
requires  significantly  more  computational  effort.  Equation  C5  can  be  used  to  find  a  relation  between  X  and 
5.  In  nondimensional  form  this  is 


-  /:21  [a(P  +  2)  +  2^]  = 


g  .  - 

The  heat  balance  integral  forms  for  eq  Cl  and  C2  are  as  follows. 


(C8) 


«l 


dTi(X,t)  dT{0,t)' 
dx  •  dx 


-U[T{ X, t) -  T (0,  7i {x, t)dx  +  7}  -^  =  0 


(C9) 


“2 


dT2{X  +  h,t)  dT2{X,t) 
dx  dx 


-[/[T2(X  +  5,t)-T2(X,t)  - 


A 

dt 


jT2ix,t)dx+T2{X  +  5)^^^^-Tf^  =  0. 


(CIO) 


These  two  energy  equations  are  summed,  along  with  eq  C3,  to  yield  an  integrated  equation  for  the  entire 
region  undergoing  temperature  changes.  The  result  is 


f.  .  r 


dt 


jpiC,  J7j(x,^)^^X  +  P2C2  J  T2{x,t)dx-p^£X  +  {p2C2-piCj)TfX-p2C2{X  +  5) 

[  {)  X 


7;+f(^+5) 


=  .  dT{0,t)  _ p^c2U[aT  +  G(X  +  5)1. 

dx 

Equation  Cl  1,  the  energy  integral  equation,  can  now  be  written  nondimensionally  as 


a 

T  =  J  K^da 

0 


= 


bi  +  b2^  -  — 
6g 

fl-— 1 
8  ) 

--^(a  +  (|))ap' 

1 

a 

fi-2l 

U  J 

+  ^21  “ 

|l  + 1  /  S'j  +  C2  J 

[(l)-i-a(P  +  l)] 

1 

(Cll) 


(C12) 


(C13) 


where 


UATi 

Gai 
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,  1  1  ^  A 

b\=-  T  +  -^  +  ^21<l> 

13 


The  derivatives  of  p  and  g  can  be  found  from  the  following  equations 


^  _  p,  _  ^5  +  aia4 

do  ^3  +  0-2^4 


^  =  g'  =  ai-a2^' 

do 


«21  ^  (t^  +  <l>)P(P  +  2) 


.22i(^|2o{|S  +  1)  +  2*] 


8  ^7 


-4=  ^  +  4li 

{  8  J 


^5  =^2i(P  +  2) 

m  =  p[^c(p  +  2]  +  2(|)|. 

The  problem  has  now  been  reduced  to  a  simple  numerical  quadrature  of  eq  C12  using  the  auxiliary  rela¬ 
tions  of  eq  C13-C15.  A  FORTRAN  program  to  carry  out  the  integration  is  listed  as  PFTSYNB.FOR  in  Ap¬ 
pendix  E. 


Phase  change  model  verification 

A  simplification  of  this  problem  can  be  solved  in  a  closed  form.  Consider  the  case  of  a  soil  initially 
thawed  at  Tf  and  with  a  zero  geothermal  gradient  G.  The  problem  is  then  one  of  a  single  phase  only  with  eq 
Cl,  Cla,b,  C3,  C4,  and  C9  governing  the  freeze  process.  The  temperature  is  chosen  as 


T  =  Tf  +  P\ 


c  „2(x-X 


X  2( 


p=-R,  /?  =  Ji+; 


The  location  of  the  freeze  interface  is  given  by 


1  K2X-^K, 
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Table  Cl.  Comparison  of  closed 
solution  (G  =  0)  and  numerical 
quadrature  (G  =  0.0001). 

St  =  0.144,  a  =  58.89  m^/yr,  = 
10°C,  U=  1  mm/yr. 


Freeze 

depth 

(m) 

Time 

(yr) 

eqC12 

Time 

(yr) 

eqC17 

Percent 

dijf'erence 

1000 

55,867 

54,778 

1.99 

2000 

206,996 

203,653 

1.64 

3000 

437,935 

428,448 

2.20 

Ky  =\  +—{\  +  RI'i) 

K2=V{1  +  Sj) 

a:3  =  (x/?(i+/?). 

Note  that  if  U  is  zero,  the  phase  change  interface  is 


2000  206,996  203,653  1.64 

x'^  = - i - 1 - -2ar.  (CIS)  3000  437,935  428,448  2.20 

l  +  R/2{l-\-R/3)  - 

This  is  identical  to  the  well  known  Stefan  solution  given  in  Lunardini  (1991).  We  may  compare  the  closed 
form  solution  (for  which  G  =  0)  with  the  numerical  quadrature  of  eq  C12  by  letting  G  be  very  small.  Table 
Cl  shows  that  the  results  are  quite  good  even  for  very  long  freeze  times. 

Method  2 

We  can  examine  the  same  problem  with  a  different  approximation  method  by  referring  to  Figure  13.  For 
region  3,  a  quasi-steady  approach  will  be  used,  leading  to  a  linear  temperature  profile.  The  basic  equations 
for  heterogenetic  growth  are  valid  except  that  the  surface  temperature  will  be  replaced  by  a  transient  func¬ 
tion  Equations  1-5  are  valid  but  the  temperature  profiles  are  changed  as  follows.  Quadratic  tempera¬ 
ture  profiles  in  regions  1  and  2  and  a  linear  temperature  in  region  3  that  satisfy  the  boundary  conditions  are 


+  {a^X--AT^ 


T2=Tf+  [g(5  +  2X)  +  2AT]^^-^-  {GX  +  AT)  ^ 


T3  =T,  +  ATiMRX  — +  1 


g 


^  =  +  /?  =  aA3(2-l/^)/( 


6[g(5  +  2X)  +  2A7’] 


M  =  -2-,  AT,  =Tf-  TXt)  =  RaATy  M  /[a,^,3  (2  -Mg)]. 

1  +  /v 

Equation  5  can  be  used  to  find  a  relation  between  X  and  8.  In  nondimensional  form  this  is 

g 

Equation  3,  the  energy  integral  equation,  can  now  be  written  nondimensionally  as 
0 

x  =  jK/Hd(^ 

0 

K  =  bi+{b2-Aa/3)^-M/{6g)-Aa/2  +  ajM^kyi{0.5/g  +  l)^[g'/g^y{2-Mg)/c  13 
+  a^{b2-Aa/3)^'  +  Mg'  /(6g^ )  -  A(p  /  3  +  0.5) . 
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(C26) 

(C27) 


The  problem  has  now  been  reduced  to  a  simple  numerical  quadrature  of  eq  C23  using  the  auxiliary  rela¬ 
tions  of  eq  C24-C26.  A  FORTRAN  program  to  carry  out  the  integration  is  listed  in  Appendix  E  as 
PFTSYN.FOR. 

This  approximation  is  inferior  to  method  1  since  the  variables  in  eq  C23  are  not  strictly  separable. 
Nevertheless,  predictions  for  modest  times  compare  quite  well  to  those  of  method  1 . 
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APPENDIX  D:  ENERGY  FLOWS  AT  THE  PERMAFROST  BASE 

The  heat  flows  at  the  base  of  the  permafrost  layer  determine  the  rate  of  movement  of  the  permafrost  bot¬ 
tom.  Examine  Figure  D1  to  clarify  the  concepts  involved.  At  any  instant  of  time,  an  amount  of  energy, 
^f(0,  is  conducted  away  from  the  phase-change  interface  through  the  frozen  layer  and  to  the  interface  from 
the  thawed  zone  2,  H  is  important  to  realize  that  the  constant  geothermal  heat  flow  =  k^G  will  only 
equal  ^u(0  at  equilibrium.  During  movement  of  the  permafrost  base,  ^u(0  can  be  greater  or  less  than  qg.  We 
can  examine  the  transient  behavior  of  these  terms.  Let 


0  gf(0  -  1  971  (X,f) 

k^G  G  dx 


Qu 


1  dT2(X,t) 
G 


In  terms  of  dimensionless  parameters 

Q{  -  7  =r 

|p[0((3  +  2)  +  2(t)]  j 

g(P  +  2)  +  2(|) 
aP 


(Dl) 

(D2) 

(D3) 


(D4) 


x  =  0 


Frozen,  1 


x  =  X(t), 
Bottom  of 
Permafrost  at  t 


Thawed,  2 


x  =  X+6, 
Depth  for 
qu  =  'lg 


Thawed,  3 

Figure  Dl.  Energy  flows  at  base  of  permafrost 


where  we  assume  that  -  ki,  i.e.,  the  thermal  conductivity  of  the  thawed  zone  is  constant.  The  results  of 
example  2  for  Prudhoe  Bay  are  tabulated  in  Table  DL  Clearly,  the  heat  flow  from  the  deep  layer  greatly  ex¬ 
ceeds  the  geothermal  heat  flow  for  much  of  the  permafrost  formation  period.  Also,  note  how  rapidly  the 
heat  flow  out  of  the  frozen  zone  Qf  drops  to  slightly  more  than  the  flow  from  the  thawed  layer  Qu* 

The  temperature  in  the  frozen  zone  adjusts  very  quickly  toward  the  equilibrium  result  of  a  linear  temper¬ 
ature  distribution;  however,  the  thermal  zone  approaches  equilibrium  very  slowly. 


Surface  temperature  increase 

Suppose  the  permafrost  has  thickness  Xq, 
is  growing  and  the  surface  temperature  ini¬ 
tially  increases  by  a  certain  amount  and  is 
held  constant  for  several  thousand  years. 
What  is  the  effect  on  the  permafrost?  The 
bottom  of  the  permafrost  will  continue  to 
grow  for  several  hundred  years  before  start¬ 
ing  to  thaw,  but  Table  Dl  assures  us  that 
this  will  be  negligible.  At  r  =  225,000  years, 
Xo  =  626.5  m  and  the  permafrost  growth 
will  only  be  27.5  cm  during  the  next  1000 
years.  Thus,  we  can  assume  that  the  perma¬ 
frost  thickness  remains  essentially  constant 
while  the  temperature  in  the  frozen  zone 


Table  Dl.  Freeze  at  Prudhoe  Bay,  Alaska. 

ATi  =  12.69°C  ({)  =  0 

k2\  =  0.5812  St  =  0.1827 

a2i  =  0.3915  Xe  ^  763.5  m 


q/  1 

Permafrost  - = - 


Time 

(yr) 

depth 

(m) 

p 

4/ 

<3/ 

“21  ,  ^ 

p(p  +  2) 

1 

4.44 

1.9060 

2.423 

159.057 

0.9251 

350 

79.9 

1.6717 

2.196 

8.986 

0.9405 

3500 

219.3 

2.3380 

1.855 

3.352 

0.9629 

35,000 

461.4 

4.8154 

1.4153 

1.6351 

0.9826 

100,000 

567.8 

7.7133 

1.2593 

1.3375 

— 

225,000 

626.5 

11.1633 

1.1792 

1.2153 

— 

775,000 

687.7 

20.3168 

1.0984 

1.1090 

— 
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adjusts  to  its  new  equilibrium  value.  Figure  D2  sketches  the  prob¬ 
lem.  Obviously,  after  infinite  time  the  new  temperature  profile  is 
as  shown.  However,  the  temperature  will  adjust  to  near  the  new 
equilibrium  in  a  relatively  short  time.  This  is  a  linear  problem  in 
non-phase-change  conduction  and  has  been  solved  by  Lachen- 
bruch  et  al.  (1982).  The  transient  temperature  is 


oo 

1  ■ 
—e  sm 

(  \ 

X 

rnt — 

(  s  sj^ 

n=\ 

n 

1  xj 

(D5) 


where 


Figure  D2.  Permafrost  equilibrium  tem¬ 
perature  profiles. 


M  = 


’ 


tc  is  a  characteristic  time  for  sensible  temperature  changes.  The  equilibrium  temperatures  are  simply 


Aq 

The  change  in  sensible  heat,  going  from  the  state  at  r  =  0  to  the  state  at  r  =  ©o,  is 

The  change  in  sensible  heat  at  any  time  t  is 

Q,=q]  {T-T„)dx  =  Q,^  +  ^(T:’ -T^X.'Z  —  l-H)"  • 

^  TC  '  n 

0 

Thus,  the  relative  change  in  sensible  heat  is 


(D6) 

(D7) 


(D8) 


(D9) 


=  (DIO) 

Qsoo  n=l  ^  ^  ^ 


Note  that  this  quantity  does  not  depend  upon  the  surface  temperatures.  The  relative  change  is  shown  in 
Table  D2.  The  sensible  heat  change  attains  93%  of  its  ultimate  value  at  t/X  =  1.0  (r  =  1666  years)  and  99%  at 
t/X  =  1.78  (t  =  2966  years).  The  sensible  heat  changes  would  be  essentially  completed  after  about  1670 
years.  From  this  time  on,  the  bottom  of  the  permafrost  would  slowly  melt. 


Change  in  frozen  zone  temperature  gradient  at  bottom  of  permafrost 

The  time  required  for  the  temperature  gradient,  in  the  frozen  zone,  to 
change  is  important  since  this  quantity  will  determine  the  rate  of  change  of  the 
permafrost  bottom  depth.  The  gradient  at  the  depth  Xq  can  be  found  from  eq 
D5  and  is 

UX  X  —  X(^  Aq  Aq 


Table  D2.  Relative  change 
in  sensible  heat. 


MX 

Qs^-Qs 

Q.oe 

0.25 

0.438 

0.50 

0.236 

1.0 

0.069 

1.50 

0.020 

1.78 

0.010 

2.0 

0.006 
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The  initial  temperature  gradient,  at  x  =  Xq,  is 


5‘'=|^(Xo,0)  =  ^^;^.  (D12) 

ax  Aq 

Note  that  this  simplified  form  of  8"  gives  virtually  the  same 
heat  flux  as  the  value  from  eq  Dl,  when  x  =  Xq-  Now 


5" -6 
5" 


rri/  rpO 

~ 


Ti-V 


l  +  2X(-l)''e-" 


(D13) 


The  time  for  a  given  change  in  the  gradient  can  be  closely  ap¬ 
proximated  as 


Table  D3.  Time  for  change  in  gradi¬ 
ent. 

7;"  =  -13.69°C;  7;' =-11.0; 


5"  -6 


5"  -6 


5 

t/X 

(years)  ^ 

1  s  ; 

0.01 

0.301 

502 

0.0314 

0.0297 

0.3421 

570 

0.0440 

0.1 

0.540 

899 

0.1022 

0.20 

1.446 

2408 

0.20 

0.212 

oo 

oo 

_ 

*  Calculated  without  eq  D14  approximation. 


t 

X 


2 


'"5"  -8^ 

5  J 

(D14) 


Table  D3  shows  the  results  for  Prudhoe  Bay. 

It  would  take  about  490  years  for  the  bottom  growth  to  cease  and  900  years  for  the  bottom  gradient  to 
change  significantly.  Thus,  the  approximations  used  in  the  derivation  of  eq  D5  are  acceptable. 


Bottom  melt 

The  temperature  in  the  frozen  zone  requires  1666  years  for  sensible  heat  adjustment,  leaving  At  =  15000  - 
1666  =  13,334  years  for  bottom  melt  (the  interglacial  is  15,000  years  long).  The  energy  balance  at  the  bottom 
of  the  permafrost  is 


X 


-Aqg. 


(D15) 


where  A  is  the  fraction  of  the  geothermal  energy  that  goes  into  melting;  it  can  exceed  1.0. 

The  heat  flow  from  the  thawed  zone  is  greater  than  the  geothermal  heat  flow,  as  we 
have  seen.  For  the  example  discussed  here,  A  =  1.179  at  the  beginning  of  thaw  and  will 
decline  towards  1  as  thaw  proceeds.  The  permafrost  thickness  Xf  after  At  years  is  given  by 


Xf-b  In 


''  Xf+b'" 


yXf^+b 


—  <3  +  Xf. 


where 


(D16) 


L  Aq^ 


q^ 


Table  D4.  Per¬ 
mafrost  thick¬ 
ness.* 


Xf 

A  (m) 

1.179  591.0 

1.0895  605.5 

LO _ 620.0 

*Final  thickness 
after  surface  tem¬ 
perature  increase, 
Prudhoe  Bay. 


The  final  permafrost  thickness  is  strongly  dependent  upon  the  value  of  A.  Table  D4  shows  values  for  the 
Prudhoe  Bay  example. 

The  heat  flow  from  the  thawed  zone  varies  continuously  during  the  thaw,  denoted  in  eq  D15  as  A{t)  qg. 
The  heat  flow  at  equilibrium  is  such  that  A  =  Ag  =  1.0.  Thus,  let  A(0  be  a  linear  function  of  X,  given  by 


A  =  (A„-Aj 


l-y, 


+  A« 


e  J 


(D16a) 
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where 


j  =  X/Xo 

Aq  =  value  at  thaw  commencement 
Aq  =  equilibrium  value 
T  —T' 

Xq  -  — - ^  =  new  equilibrium  permafrost  thickness. 

/CyfG 

Equation  D15  then  has  the  following  solution 


In 


^  «1  +  «2  ”  1  J 


^  yj  +^4  ^  ^1+^5 

>'f  +<^5  JU  +  ^4  y 


+  —  0 


(D17) 


where 

_  (-^0  ^e) 


^2  3^e 


K-4) 

(l-3'e) 


(24  = 


^2  ~'^3 

2a, 


«5  = 


+^3 

20, 


tf  = 


kf[Tf-T:)At 

LXl 


and 


is  the  permafrost  thickness  after  Dr  years  (Table  D5). 

Note  that  these  results  agree  quite  well  with  the  values  with  a  constant  thawed  zone  heat  flow,  i.e.,  con¬ 
stant  value  of  A(r).  For  this  case,  15,000  years  is  nearly  enough  time  to  thaw  back  to  the  new  equilibrium 
thickness  of  601.5  m. 


Table  D5.  Permafrost 
thickness  after  thaw. 


A,,=\. 

1792. 

Ae 

y/ 

X/fmj 

1.0 

0.9678 

606.4 

1.170 

0.945 

592.0 

1.179 

0.9434 

591.0 
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APPENDIX  E:  FORTRAN  PROGRAMS  FOR  NUMERICAL  QUADRATURE 

OF  ENERGY  EQUATION 


Program  PERM 

C  $DECLARE 

C  $DEBUG 

CCC  10  REM  'k-k'k-k'k'k'k'k'k'k  PERM  'k'k'k'k-k-k'k'k:k'k 
CCC  REAL  A21 , BET , C21 , EP , K , K21 , PHI , R21 , R , ST , Y , B1 , B2 , B3 , B4 
CCC  REAL  D , DELS , DELTATO , DELTATl , G , GG , S IGI , S IGF , XI 
CCC  REAL  Al.A2,A3,A4,A5,TAU,X,T,B,BPP,GP,GPP,F,H,FP 
IMPLICIT  DOUBLE  PRECISION  (A-H , K-M, 0-Z) 

COMMON  /DATA/  A21 , BET . K , K21 , PHI , R21 , R, ST , Y , B1 , B2 , B3 
character  PRNTR*12 
integer  i , j , n 
c 

PRNTR  =  'OUT. DAT' 

2  OPEN(9,FILE-PRNTR,  STATUS =' NEW' ) 

c 

C  INPUT  VALUES  ******-A"**^»r'A-**^*^t-5^* 

ST  =  0.143960 
EP  =  0.379 
KI  =  5.45 
KW  =  1.34 
KG  =  10.360 
GAM  =  .9825 

C  CALCULATED  VALUES  *********>Sr* 

KWI  -  KW/KI 

K21  =  ((KWI)**EP)*((GAM)*^(1.0-EP)) 

C21  =  (.9357  +  1.243*EP)/(.9155+.0845*EP) 

A21  ==  K21/C21 

R21  -  (2.6-1.6*EP)/(2.6-1.69*EP) 

Kl  =  ((KI)^*EP)*((KG)**(1.0-EP)) 

Cl  =  .4202+.0388*EP 
ALl  -  3.1536*K1/C1 
L  =  79.71*EP 

DELTATO  =0.0 
DELTATl  =  L*ST/C1 

C  INCREMENT  FOR  SIMPSON  ***'A-*** 


DELS  =  0.00005 
D  =  DELS/100. 

GG  =  0.0286 

BET  =1.0 

R  =  0.33330 

PHI  =  DELTATO/DELTATl 

SIGI  =  1.0/K21 

SIGF  =  0.9*SIGI 

XI  =  DELTATl/ (GG*K21) 

B3  =  K21/A21 

B1  =  -(1. 0/3.0  4-  l./ST  +  B3*PHI) 
B2  =  -B3*PHI/3.0 
B4  =  SIGI/DELS 
N  =  B4 

130  Write(9,3000) 
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3000  formatClx/  TIME  OF  PERMAFROST  FORMATION  ') 


137  Write(9,3015)  EP 

3015  forraat(lx,'  VOID  RATIO  EPSILON  =  ',F9.4) 

140  Write(9,3001)  ST 

3001  format(lx/  STEFAN  NUMBER,  C1*DELTAT1/L  =  \F9.4) 

141  Write(9,3004)  R 

3004  format(lx,  '  SUSPECT  CONSTANT  =  \F9.4) 

150  Write(9,3002)  K21 

3002  format(lx,'  THAW/FROZEN  CONDUCTIVITY  RATIO  =  ',F9.4) 

138  Write(9,3016)  C21 

3016  format(lx,'  THAW/FROZEN  HEAT  CAPACITY  RATIO  -  ',F9.4) 

160  Write(9,3003)  A21 

3003  format(lx,'  THAW/FROZEN  DIFFUSIVITY  RATIO  =  ',F9.4) 

Write(9,3013)  R21 

3013  format(lx/  THAW/FROZEN  DENSITY  RATIO  «  ',F9.4) 

Writ:e(9,5)  KG 

5  formatClx, '  SOIL  SOLIDS  KG  MCAL/(S -CM-C)  =  ',F9.4) 
Write(9,3021)  K1 

3021  format(lx.'  FROZEN  K  MCAL/(S- CM-C)  =  ',F9.4) 

Write(9,3022)  Cl 

3022  formatdx,'  FROZEN  C  CAL/(CM**3 -C)  =  ',F9.4) 

Write(9,3023)  ALl 

3023  format(lx,^  FROZEN  ALPHA  M*’*: 2 /YEAR  =  dF9.4) 

Write(9,3024)  L 

3024  format(lx.'  LATENT  HEAT  CAL/CC  =  ',F9.4) 

Write (9. 3014)  DELTATO 

3014  format(lx,'  (TO  -  TF)  DEGREE  C  =  ',F9.4) 

180  Write(9,3005)  DELTATl 

3005  format(lx/  (TF  -  TS)  DEGREE  C  =  ',F9.4) 

191  Write(9,3007)  GG 

3007  format(lx,'  GEOTHERMAL  GRADIENT  DEG  C/M  «  ',F9.4) 

200  Urite(9,3008)  D 

3008  format(lx,'  INCREMENT  FOR  SIMPSON  “  ',F9.4) 

201  Write(9.3009)  PHI 

3009  format(lx,'  PHI  =  (TO  -  TF)/(TF  -  TS)  =  ',F9.4) 

202  Write(9,3010)  SIGI 

3010  format(lx,'  EQUIBRIUM  FREEZE  DEPTH  =  ',F9.4) 

203  Write(9,3019)  SIGF 

3019  format(lx,'  90  %  EQUIBRIUM  FREEZE  DEPTH  =  ',F9.4) 
210  Write(9.3011)  XI 

3011  format(lx,'  EQUIBRIUM  FREEZE  DEPTH  METERS  =  ',F9.4) 


TAU  =  0.0 
X  -  0,0 

222  Write(9,*) 

252  Write(9,*)  '  TIME  SIGMA  BETA 

253  Write(9,2001)  TAU,  X.  BET 
2001  format (lx , F13 . 4 , 2F10 . 4) 

260  DO  80,  I-1,N 
TS  -  X 

TE  =  X  +  DELS 

IF(X  .EQ.  0.0)  GO  TO  10 

Y  =  TS 
CALL  BETA 
KS  K 

11  Y  “  TS  +D 
SO  =  0.0 
DO  20  J  =1,50 
CALL  BETA 
SO  =  SO  +  K 
Y=Y+2*D 
20  CONTINUE 

SOO  =  4*S0 

Y  =  TS-I-2.0*D 
SE  =  0.0 

DO  30  JJ  =  1,49 
CALL  BETA 
SE  =  SE+K 

Y  =  Y+2.0*D 

30  CONTINUE 

SOE  =  2.0*SE 

Y  =  TE 


40 


CALL  BETA 
KE  =  K 
SI  =  KS+KE 

T  -  D*(Sl+S00+S0E)/3.0 
TAU  =  TAU+T 
X  -  X  +  DELS 
131  Write(9,*) 

133  Write(9,2001)  TAU,  X,  BET 

80  CONTINUE 

GO  TO  1850 

10  KS  =  0.0 

GO  TO  11 
1850  END 

SUBROUTINE  BETA 

IMPLICIT  DOUBLE  PRECISION  (A-H , K-M , 0- Z) 

COMMON  /DATA/  A21 , BET , K, K21 , PHI , R21 , R , ST , Y , Bl , B2 , B3 

40  BET  -  BET 

M  =  Y’^(BET+2.0)+2,0*PHI 

G  =  A21*(Y+PHI)/(BET*M)+1.0 

GP  =  -A21*(Y+PHI)*(M+BET*Y)/((BET*M)**2) 

B  =  BET/G  -  M*K21  -  2 , 0*R21*BET’*^(G- 1 . 0)/ST 

BP  =  (1.0-BET*GP/G)/G  -  K21*Y  +2 .0*R21*A21*Y*(Y+PHI)/(ST*(M**2) ) 
BETl  =  BET  -  B/BP 

IF(ABS(BET1-BET)  .LT.  0.01)  GO  TO  50 
BET  =  BETl 
GO  TO  40 
50  BET  -  BETl 

M  =  Y*(BET+2.0)+2.0*PHI 

G  -  A21*(Y+PHI)/(BET*M)+1.0 

Al  -  A21*(1.0-(Y+PHI)*(BET+2.0)/M)/(M*BET) 

A2  -  2.0*A21*(Y+PHI)*(Y*(BET+1.0)+PHI)/((M*BET)**2) 

A3  -  1.0/G  -  2.0*R21*(G-1.0)/ST  -  K21*Y 
A4  =  (2.0*R21/ST+1.0/(G*^2))*BET 
A5  =  K21*(BET+2.0) 

BPP  =  (A5+Al*A4)/(A3+A2*A4) 

GPP  =  Al  -  A2*BPP 

F  =  Bl+(B2-B3*Y/3)*BET-1.0/(6*G)-B3*Y/2.0 

H  =  (1.0/G-2.0)/Y+K21 

FP  -  (B2-B3*Y/3)*BPP+GPP/(6*(G**2))-B3*(BET/3+.5) 

K  =  (F4-Y*FP)/H 

RETURN 

END 


Program  PFTSYN 

C  $DECLARE 

C  $DEBUG 

Ccc  10  REM  PFTSYN  i('k'k-:k'k-k-k-k-k-k 

CCC  REAL  A21 , BET , C21 , EP , K . K21 , K13 , PHI . R21 , R , ST . Y , Bl , B2 , B3 , B4 

CCC  REAL  D . DELS , DELTATO . DELTATl ,G,GG,PSI,SIGD,SIGI,SIGF,SIGL,U.XI 

CCC  REAL  A,Al,A2,A3,A4,A5,TAU,X.XD,XP,XT,T,B,BPP.BPA,GP,GPP,F,H.FP 

CCC  REAL  S1,S2,S3,MN,H1,Q1,Q2,Q3,Q4,Q5,Q6,Q12,Q13,Q14.Q15 

IMPLICIT  DOUBLE  PRECISION  (A-H , K-M, 0-Z) 

COMMON  /DATA/  A , A21 , BET , K , K21 , K13 , XD , PHI , PSI , R21 , R , ST , Y , Bl . B2 . B3 
character  PRNTR*12 
integer  i , j ,n 
c 

PRNTR  =  'OUT. DAT ^ 

2  OPEN(9.,FILE-PRNTR,  STATUS^ ' NEW' ) 

c 

C  'k'k'k'k'k'kii'k  INPUT  VALUES  'k'k'k'k'k'k^k-k'k'k-kik-k-k-k-k-k'k 
ST  -  .143960250 
EP  =  0.379 
KI  =  5.45 
KU  =  1.34 
KG  =  10.360 
GAM  =  .9825 
U  =  10.0 
K13  =  1.0 
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c  CALCULATED  VALUES  *********^* 

KWI  =  KW/KI 

K21  =  ((KWI)**EP)*((GAM)*^(1.0-EP)) 

C21  =  (.9357  +  1.243*EP)/(.9155+.0845*EP) 

A21  =  K21/C21 

R21  =  (2.6-1.6*EP)/(2.6-1.69*EP) 

K1  =  ((KI)**EP)*((KG)’^*(1,0-EP)) 

Cl  =  .4202+.0388*EP 
ALl  =  3.1536*K1/C1 
L  -  79.71*EP 

C  K21  -  1.0 

C  C21  =  1.0 

C  A21  =  1.0 

C  R21  =  1.0 

C  ST  =  1000. 

A«C21 

DELTATO  =0.0 
DELTATl  =  L*ST/C1 
C  DELTATl  =10.0 

C  -k-kic^-k-k  INCREMENT  FOR  SIMPSON  kkkkkkk 

DELS  =  0.00005 
D  =  DELS/100. 

GG  =  0.0286 

PSI  =  U’t.001*DELTATl/(ALl*GG) 

BET  =1.0 

PHI  =  DELTATO/DELTATl 
SIGI  =  1.0/K21 
SIGF  =  0.9*SIGI 
SIGL  =  SIGI-0.01 
XI  =  DELTATl/ (GG*K21) 

B3  =  K21/A21 

B1  =  -(1. 0/3.0  +  l./ST  +  B3*PHI) 

B2  =  -B3*PHI/3.0 
B4  =  SIGI/DELS 
N  =B4 

130  Write(9,3000) 

3000  format(lX,'  TIME  OF  SYNGENETIC  PERMAFROST  FORMATION  ') 


137  Write(9,3015)  EP 

3015  format(lx,'  VOID  RATIO  EPSILON  =  '.F9.4) 

140  Write(9,3001)  ST 

3001  format(lx,'  STEFAN  NUMBER,  C1*DELTAT1/L  =  ',F9.4) 

141  Write(9 , 3004)  U 

3004  formatClx,'  DEPOSITION  RATE  MM/YR  =  ',F9.4) 

129  Write(9,2999)  PSI 

2999  format(lx,'  DEPOSITION  PARAMETER  PSI  =  ',F9.6) 

150  Write(9,3002)  K21 

3002  format(lx,'  THAW/FROZEN  CONDUCTIVITY  RATIO  =  ^F9.4) 

138  Write(9.3016)  C21 

3016  format(lx,'  THAW/FROZEN  HEAT  CAPACITY  RATIO  =  ',F9.4) 

160  Write(9,3003)  A21 

3003  format(lx,'  THAW/FROZEN  DIFFUSIVITY  RATIO  =  ',F9.4) 

Write(9,3013)  R21 

3013  format(lx,'  THAW/FROZEN  DENSITY  RATIO  =  ',F9.4) 

Write(9.5)  KG 

5  format(lx,'  SOIL  SOLIDS  KG  MCAL/(S -CM-C)  =  ',F9.4) 

Write(9,3021)  Kl 

3021  format(lx/  FROZEN  K  MCAL/(S- CM-C)  =  ',F9.4) 

WrLte(9,3022)  Cl 

3022  formatdx,'  FROZEN  C  CAL/(CM**3 - C)  =  ',F9.4) 

Write(9,3023)  ALl 

3023  format(lx,^  FROZEN  ALPHA  M**2/Y£AR  =  ^F9.4) 

Write(9,3024)  L 

3024  format(lx/  LATENT  HEAT  CAL/CC  =  '.F9.4) 

Write (9, 3014)  DELTATO 

3014  format(lx.'  (TO  -  TF)  DEGREE  C  =  ',F9.4) 

180  Write(9,3005)  DELTATl 

3005  format(lx,^  (TF  -  TS)  DEGREE  C  =  ',F9.4) 

191  Write(9,3007)  GG 

3007  format(lx,'  GEOTHERMAL  GRADIENT  DEG  C/M  =  ',F9.4) 

200  Write(9,3008)  DELS 

3008  forniat(lx,'  INCREMENT  FOR  SIMPSON  =  \F9.6) 

201  Write(9, 3009)  PHI 
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3009  formatClx/  PHI  =  (TO  -  TF)/(TF  -  TS) 

202  Write(9,3010)  SIGI 

3010  format(lx,'  EQUIBRIUM  FREEZE  DEPTH 

203  Write(9,3019)  SIGF 

3019  format(lx/  90  %  EQUIBRIUM  FREEZE  DEPTH 
210  Write(9.30Il)  XI 

3011  format(lx,'  EQUIBRIUM  FREEZE  DEPTH  METERS 
TAU  =0.0 

X  =  0.0 
XD  =  0.0 
XP  =  0.01 


' ,F9.4) 
^ ,F9.4) 
' ,F9.4) 
' ,F9.4) 


XT-X+XD 

222  Write(9,*) 

252  Write(9,*)  '  TIME  SIGMA  BETA  SIGD 

253  Write(9,2001)  TAU,  X.  BET,XD,XT 
2001  format (lx , F13 . 4 , 4F10 . 4) 


260  DO  80  1=1, N 
TS  =  X 

TE  =  X  +  DELS 

IF(X  .EQ.  0.0)  GO  TO  10 

Y  =  TS 
CALL  BETA 
KS  =  K 

11  Y  =  TS  +D 
SO  =  0.0 
DO  20  J  =1,50 
CALL  BETA 

50  =  SO  +  K 
Y=Y+2’^D 

20  CONTINUE 

SOO  =  4*S0 

Y  =  TS+2.0*D 
SE  =  0.0 

DO  30  JJ  =  1,49 
CALL  BETA 
SE  =  SE+K 

Y  =  Y+2.0*D 

30  CONTINUE 

SOE  =  2.0*SE 

Y  =  TE 
CALL  BETA 
KE  =  R 

51  =  KS+KE 

T  =  D*(S1+SOO+SOE)/3.0 
IF(T  .LT.  0.0)  GO  TO  131 
TAU  =  TAU+T 
X  =  X  +  DELS 
XD  =  XD+PSI*T 
XT  =  XD+X 

IF(XT  .GT.  SIGL)  GO  TO  131 
IF(X  .LT.  XP)  GO  TO  80 
XP  =XP  +0.01 
131  Vrite(9,*) 

133  Write(9,2001)  TAU,  X,  BET,  XD,  XT 

IF(XT  .GT.  SIGL)  GO  TO  1850 
IF(T  .LT.  0.0)  GO  TO  1850 
80  CONTINUE 

GO  TO  1850 
10  KS  =  0.0 

GO  TO  11 
1850  END 


SIGT 


SUBROUTINE  BETA 

IMPLICIT  DOUBLE  PRECISION  (A-H , K-M , 0- Z) 

COMMON  /DATA/  A , A21 , BET , K , K21 , K13 , XD , PHI , PS  I , R21 , R , ST , Y , B1 , B2 . B3 

40  BET  =  BET 

M  =  Y*(BET+2.0)+2.0*PHI 
G  =  A21*(Y+PHI)/(BET*M)+1.0 
R  =  K13*XD*(2. -l./G)/Y 
GP  =  -A21*(Y+PHI)*(M+BET*Y)/((BET*M)**2) 

B  =  BET/G  -  M*K21  -  2 . 0*R21*BET*(G - 1 . 0)/ST-R*BET/(  ( 1 . +R)  ^"G) 

PI  =  1.0/((1.+R)^(2.*G-1.)) 

P  =  R/((1.0+R)*G) 


43 


BPA-P+P*BET*GP* ( PI - 1 . ) /G 

BP=(1 . 0-BET*GP/G) /G-K21*Y+2 . 0*R21*A21*Y*(Y+PHI)/(ST*(M**2) ) -BPA 
BETl  =  BET  -  B/BP 

IF(ABS(BET1-BET)  .LT.  0.01)  GO  TO  50 
BET  =  BETl 
GO  TO  40 
50  BET  -  BETl 

M  -  Y*(BET+2.0)+2.0*PHI 
G  =  A21*(Y+PHI)/(BET*M)+1.0 
R  =  K13*XD*(2. -l./G)/Y 

Al  -  A21*(1.0-(Y+PHI)*(BET+2.0)/M)/(M*BET) 

A2  =  2.0*A21*(Y+PH1)*(Y*(BET+1.0)+PHI)/((M*BET)**2) 
Sl=-A*(BET*(Y+PHI)/3 . +Y/2 . +PHI ) - 1 . /ST 

52- K13*(2.-l./G)/Y 

53- 0,5/G+1.0 
MN=l./(l-0+R) 

Hl=MN*(l./G-2. )/Y+K21-Y*S3*(MN**2)*PSI*S2/3. 

Q12-S  3* (MN**2 ) *XD*K13/ ( 3 . *G ) 

Ql=Sl-MN*S3/3.+<!12*(l.  -2.*G)A-Y*A*(BET/3.  +  .5) 

Q2=Q1 2 /G+MN*Y/ ( 6 . * (G**2 ) ) 

Q3~Y*A*(Y+PHI)/3. 

Q13=PSI*S2/H1 

Q14=BET*(MN**2)/G 

q4=.q14*(Q13*Q1+XD*K13*(1 . -2 . *G)/( (Y**2)*G) ) 
q5=.-Q14*(Q13*Q2+XD*K13/(Y*(G**2))) 
q6=Q14*Q3*Q13 
Q15=2 .*R21/ST 

A3  =  MN/G  -  Q15*(G-1.0)  -  K21*Y+Q6 
A4  -  Q15*BET+1.0/(G**2)*BET*MN-Q5 
A5  =  K21*(BET+2.0)-Q4 
BPP  =  (A5+A1*A4)/(A3+A2*A4) 

GPP  =  Al  -  A2*BPP 

K  -  (Q1+Q2*GPP-Q3*BPP)/H1 

RETURN 

END 
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